# Directory and Data Import
setwd("/home/joshk/git_repositories/BrookTrout_TG/Data")
BTChallenge <- read.csv("BT Challenge_DegreeDay.csv")
{
BTChallenge$Offspring.Acc <- as.factor(BTChallenge$Offspring.Acc)
BTChallenge$Mother.Acc <- as.factor(BTChallenge$Mother.Acc)
BTChallenge$Father.Acc <- as.factor(BTChallenge$Father.Acc)
BTChallenge$Mother.ID <- as.factor(BTChallenge$Mother.ID)
BTChallenge$Father.ID <- as.factor(BTChallenge$Father.ID)
BTChallenge$Subject.ID <- as.factor(BTChallenge$Subject.ID)
BTChallenge$Family <- as.factor(BTChallenge$Family)
BTChallenge$Parent.Group <- as.factor(BTChallenge$Parent.Group)
}
options(na.action = na.omit)
## Loading in necessary packages
library("easypackages")
libraries(
"broom", "car", "dotwhisker", "dplyr", "ggplot2", "grid",
"gridExtra", "gtable", "itsadug", "knitr", "lme4", "lmerTest",
"mgcv", "MuMIn", "nlme", "purrr", "rmarkdown", "splines"
)
# Adding font
library("showtext")
font_add_google(name = "Noto Sans", family = "Noto Sans")
## Additional functions
strip.lme <- function(model.name) {
coef.names <- c(names(model.name$coefficients$fixed))
tTab.data <- data.frame(summary(model.name)$tTable[, c(1, 2, 4)], row.names = NULL)
tTab.data <- cbind(coef.names, tTab.data)
return(tTab.data)
}
sig.lme <- function(model.name, sig.level = 0.1) {
sig.data <- as.data.frame(summary(model.name)$tTable[c(which(as.numeric(paste(summary(model.name)$tTable[, 5])) <= sig.level)), ])
return(sig.data)
}
Assessing parameter distributions by plots
Unique.Mass <- distinct(BTChallenge, Mass, .keep_all = TRUE) %>%
select(Subject.ID, Mass, Offspring.Acc, Mother.Acc, Father.Acc, Family)
ggplot(data = Unique.Mass, aes(x = Mass)) +
geom_histogram(colour = "black", fill = "dodgerblue2", bins = 20) +
theme_bw() +
my.theme.diag +
geom_vline(
xintercept = with(na.omit(Unique.Mass), mean(Mass)), size = 1.25,
colour = "black", linetype = "dotted"
) +
xlab("Mass (g)") +
geom_density(alpha = 0.5, aes(y = 1 / 5 * ..count..), fill = "dodgerblue2")
ggplot(data = Unique.Mass, aes(
x = Mass, fill = as.factor(Offspring.Acc),
colour = as.factor(Offspring.Acc)
)) +
geom_histogram(bins = 20) +
theme_bw() +
my.theme.diag +
geom_vline(
xintercept = with(na.omit(Unique.Mass), mean(Mass)), size = 1.25,
colour = "black", linetype = "dotted"
) +
xlab("Mass (g)") +
geom_density(alpha = 0.5, aes(y = 1 / 4.2 * ..count..)) +
scale_fill_manual(
values =
c("dodgerblue3", "maroon1")
) +
scale_colour_manual(
values =
c("black", "black"), guide = FALSE
) +
guides(fill = guide_legend(
title =
"Offspring Acclimation \n Temperature"
))
## Strong overlap
ggplot(data = Unique.Mass, aes(
x = Mass, fill = as.factor(Mother.Acc),
colour = as.factor(Mother.Acc)
)) +
geom_histogram(bins = 20) +
theme_bw() +
my.theme.diag +
geom_vline(
xintercept = with(na.omit(Unique.Mass), mean(Mass)), size = 1.25,
colour = "black", linetype = "dotted"
) +
xlab("Mass (g)") +
geom_density(alpha = 0.5, aes(y = 2 / 7 * ..count..)) +
scale_fill_manual(
values =
c("slateblue", "lightpink")
) +
scale_colour_manual(
values =
c("black", "black"), guide = FALSE
) +
guides(fill = guide_legend(
title =
"Maternal Acclimation \n Temperature"
))
ggplot(data = Unique.Mass, aes(
x = Mass, fill = as.factor(Father.Acc),
colour = as.factor(Father.Acc)
)) +
geom_histogram(bins = 20) +
theme_bw() +
my.theme.diag +
geom_vline(
xintercept = with(na.omit(Unique.Mass), mean(Mass)), size = 1.25,
colour = "black", linetype = "dotted"
) +
xlab("Mass (g)") +
geom_density(alpha = 0.5, aes(y = 2 / 7 * ..count..)) +
scale_fill_manual(
values =
c("seagreen3", "magenta3")
) +
scale_colour_manual(
values =
c("black", "black"), guide = FALSE
) +
guides(fill = guide_legend(
title =
"Paternal Acclimation \n Temperature"
))
## No obvious distinction between offspring from warm and cold mothers, or fathers.
Exploring metabolic responses to thermal challenge across groupings. Error bars represent +/- 95% CI.
Challenge.Sum <- BTChallenge %>%
group_by(Challenge.Temp) %>%
summarize(
Mean = mean(MO2.NMS, na.rm = TRUE), Upper.CI = mean(MO2.NMS,
na.rm = TRUE
) + 1.96 * (sd(MO2.NMS, na.rm = TRUE) / sqrt(length(na.omit(MO2.NMS)))),
Lower.CI = mean(MO2.NMS, na.rm = TRUE) - 1.96 * (sd(MO2.NMS, na.rm = TRUE) /
sqrt(length(na.omit(MO2.NMS))))
)
## All offspring grouped.
ggplot(Challenge.Sum, aes(x = Challenge.Temp, y = Mean)) +
theme_bw() +
my.theme +
geom_errorbar(
mapping = aes(
x = Challenge.Temp,
ymin = c(Lower.CI), ymax = c(Upper.CI)
), width = 0.2, size = 1,
color = "dodgerblue3"
) +
geom_point(size = 1, colour = "maroon1") +
xlab("Challenge Temperature \n (°C)") +
ylab("Metabolic Rate (mg O2/h)")
# Robust elevation at 25 °C.
Challenge.Sum.Off <- BTChallenge %>%
group_by(Challenge.Temp, Offspring.Acc) %>%
summarize(
Mean = mean(MO2.NMS, na.rm = TRUE), Upper.CI = mean(MO2.NMS,
na.rm = TRUE
) + 1.96 * (sd(MO2.NMS, na.rm = TRUE) / sqrt(length(na.omit(MO2.NMS)))),
Lower.CI = mean(MO2.NMS, na.rm = TRUE) - 1.96 * (sd(MO2.NMS, na.rm = TRUE) /
sqrt(length(na.omit(MO2.NMS))))
)
## Split by offspring acclimation temperature.
ggplot(as.data.frame(Challenge.Sum.Off), aes(
x = Challenge.Temp,
y = Mean, group = Offspring.Acc, colour = Offspring.Acc
)) +
theme_bw() +
my.theme +
geom_errorbar(
mapping = aes(
x = Challenge.Temp,
ymin = c(Lower.CI), ymax = c(Upper.CI)
), width = 0.2, size = 1,
position = position_dodge(width = 0.5)
) +
geom_point(size = 1, colour = "black", position = position_dodge(width = 0.5)) +
xlab("Challenge Temperature \n (°C)") +
ylab("Metabolic Rate (mg O2/h)") +
theme(legend.position = "top") +
scale_colour_manual(values = c("dodgerblue", "maroon1")) +
guides(colour = guide_legend(title = "Offspring Acclimation \n Temperature"))
# Some divergence, but does not appear systematic.
Challenge.Sum.Mat <- BTChallenge %>%
group_by(Challenge.Temp, Mother.Acc) %>%
summarize(
Mean = mean(MO2.NMS, na.rm = TRUE), Upper.CI = mean(MO2.NMS,
na.rm = TRUE
) + 1.96 * (sd(MO2.NMS, na.rm = TRUE) / sqrt(length(na.omit(MO2.NMS)))),
Lower.CI = mean(MO2.NMS, na.rm = TRUE) - 1.96 * (sd(MO2.NMS, na.rm = TRUE) /
sqrt(length(na.omit(MO2.NMS))))
)
# Split by maternal acclimation temperature
ggplot(as.data.frame(Challenge.Sum.Mat), aes(
x = Challenge.Temp,
y = Mean, group = Mother.Acc, colour = Mother.Acc
)) +
theme_bw() +
my.theme +
geom_errorbar(
mapping = aes(
x = Challenge.Temp,
ymin = c(Lower.CI), ymax = c(Upper.CI)
), width = 0.2, size = 1,
position = position_dodge(width = 0.5)
) +
geom_point(size = 1, colour = "black", position = position_dodge(width = 0.5)) +
xlab("Challenge Temperature \n (°C)") +
ylab("Metabolic Rate (mg O2/h)") +
theme(legend.position = "top") +
scale_colour_manual(values = c("slateblue", "lightpink")) +
guides(colour = guide_legend(title = "Maternal Acclimation \n Temperature"))
# Perhaps a very subtle divergence.
Challenge.Sum.Pat <- BTChallenge %>%
group_by(Challenge.Temp, Father.Acc) %>%
summarize(
Mean = mean(MO2.NMS, na.rm = TRUE), Upper.CI = mean(MO2.NMS,
na.rm = TRUE
) + 1.96 * (sd(MO2.NMS, na.rm = TRUE) / sqrt(length(na.omit(MO2.NMS)))),
Lower.CI = mean(MO2.NMS, na.rm = TRUE) - 1.96 * (sd(MO2.NMS, na.rm = TRUE) /
sqrt(length(na.omit(MO2.NMS))))
)
# Split by Paternal Acclimation Temperature.
ggplot(as.data.frame(Challenge.Sum.Pat), aes(
x = Challenge.Temp,
y = Mean, group = Father.Acc, colour = Father.Acc
)) +
theme_bw() +
my.theme +
geom_errorbar(
mapping = aes(
x = Challenge.Temp,
ymin = c(Lower.CI), ymax = c(Upper.CI)
), width = 0.2, size = 1,
position = position_dodge(width = 0.5)
) +
geom_point(size = 1, colour = "black", position = position_dodge(width = 0.5)) +
xlab("Challenge Temperature \n (°C)") +
ylab("Metabolic Rate (mg O2/h)") +
theme(legend.position = "top") +
scale_colour_manual(values = c("seagreen3", "magenta3")) +
guides(colour = guide_legend(title = "Paternal Acclimation \n Temperature"))
# Subtle effect at upper temperatures.
Challenge.Sum.OffMat <- BTChallenge %>%
group_by(Challenge.Temp, Offspring.Acc, Mother.Acc) %>%
summarize(
Mean = mean(MO2.NMS, na.rm = TRUE), Upper.CI = mean(MO2.NMS,
na.rm = TRUE
) + 1.96 * (sd(MO2.NMS, na.rm = TRUE) / sqrt(length(na.omit(MO2.NMS)))),
Lower.CI = mean(MO2.NMS, na.rm = TRUE) - 1.96 * (sd(MO2.NMS, na.rm = TRUE) /
sqrt(length(na.omit(MO2.NMS))))
)
OffMat.Group <- c()
for (i in 1:nrow(Challenge.Sum.OffMat)) {
if (Challenge.Sum.OffMat$Offspring.Acc[i] == 1 & Challenge.Sum.OffMat$Mother.Acc[i] == 1) {
OffMat.Group[i] <- ("1")
} else if (
Challenge.Sum.OffMat$Offspring.Acc[i] == 1 & Challenge.Sum.OffMat$Mother.Acc[i] == 2) {
OffMat.Group[i] <- ("2")
} else if (
Challenge.Sum.OffMat$Offspring.Acc[i] == 2 & Challenge.Sum.OffMat$Mother.Acc[i] == 1) {
OffMat.Group[i] <- ("3")
} else if (
Challenge.Sum.OffMat$Offspring.Acc[i] == 2 & Challenge.Sum.OffMat$Mother.Acc[i] == 2) {
OffMat.Group[i] <- ("4")
}
}
Challenge.Sum.OffMat$OffMat <- OffMat.Group
# Parsing by Maternal and Offspring Acclimation Temperature
ggplot(as.data.frame(Challenge.Sum.OffMat), aes(
x = Challenge.Temp,
y = Mean, group = OffMat, shape = Mother.Acc, colour = Offspring.Acc
)) +
theme_bw() +
my.theme +
geom_errorbar(
mapping = aes(
x = Challenge.Temp,
ymin = c(Lower.CI), ymax = c(Upper.CI)
), width = 0.2, size = 1,
position = position_dodge(width = 0.5)
) +
geom_point(size = 1.5, colour = "black", position = position_dodge(width = 0.5)) +
xlab("Challenge Temperature \n (°C)") +
ylab("Metabolic Rate (mg O2/h)") +
theme(legend.position = "top") +
scale_colour_manual(values = c("black", "gold2")) +
guides(colour = guide_legend(title = "Offspring Acclimation \n Temperature")) +
guides(shape = guide_legend(title = "Maternal Acclimation \n Temperature")) +
theme(legend.title = element_text(size = 10))
Challenge.Sum.OffPat <- BTChallenge %>%
group_by(Challenge.Temp, Offspring.Acc, Father.Acc) %>%
summarize(
Mean = mean(MO2.NMS, na.rm = TRUE), Upper.CI = mean(MO2.NMS,
na.rm = TRUE
) + 1.96 * (sd(MO2.NMS, na.rm = TRUE) / sqrt(length(na.omit(MO2.NMS)))),
Lower.CI = mean(MO2.NMS, na.rm = TRUE) - 1.96 * (sd(MO2.NMS, na.rm = TRUE) /
sqrt(length(na.omit(MO2.NMS))))
)
OffPat.Group <- c()
for (i in 1:nrow(Challenge.Sum.OffPat)) {
if (Challenge.Sum.OffPat$Offspring.Acc[i] == 1 & Challenge.Sum.OffPat$Father.Acc[i] == 1) {
OffPat.Group[i] <- ("1")
} else if (
Challenge.Sum.OffPat$Offspring.Acc[i] == 1 & Challenge.Sum.OffPat$Father.Acc[i] == 2) {
OffPat.Group[i] <- ("2")
} else if (
Challenge.Sum.OffPat$Offspring.Acc[i] == 2 & Challenge.Sum.OffPat$Father.Acc[i] == 1) {
OffPat.Group[i] <- ("3")
} else if (
Challenge.Sum.OffPat$Offspring.Acc[i] == 2 & Challenge.Sum.OffPat$Father.Acc[i] == 2) {
OffPat.Group[i] <- ("4")
}
}
# Now parsing by paternal and offspring acclimation temperatures
Challenge.Sum.OffPat$OffPat <- OffPat.Group
ggplot(as.data.frame(Challenge.Sum.OffPat), aes(
x = Challenge.Temp,
y = Mean, group = OffPat, shape = Father.Acc, colour = Offspring.Acc
)) +
theme_bw() +
my.theme +
geom_errorbar(
mapping = aes(
x = Challenge.Temp,
ymin = c(Lower.CI), ymax = c(Upper.CI)
), width = 0.2, size = 1,
position = position_dodge(width = 0.5)
) +
geom_point(size = 1.5, colour = "black", position = position_dodge(width = 0.5)) +
xlab("Challenge Temperature \n (°C)") +
ylab("Metabolic Rate (mg O2/h)") +
theme(legend.position = "top") +
scale_colour_manual(values = c("dodgerblue3", "maroon1")) +
guides(colour = guide_legend(title = "Offspring Acclimation \n Temperature")) +
guides(shape = guide_legend(title = "Paternal Acclimation \n Temperature")) +
theme(legend.title = element_text(size = 10))
# Possible, subtle interaction.
Producing three-way “interaction” plots
Challenge.Sum.All <- BTChallenge %>%
group_by(Challenge.Temp, Offspring.Acc, Father.Acc, Mother.Acc) %>%
summarize(
Mean = mean(MO2.NMS, na.rm = TRUE), Upper.CI = mean(MO2.NMS,
na.rm = TRUE
) + 1.96 * (sd(MO2.NMS, na.rm = TRUE) / sqrt(length(na.omit(MO2.NMS)))),
Lower.CI = mean(MO2.NMS, na.rm = TRUE) - 1.96 * (sd(MO2.NMS, na.rm = TRUE) /
sqrt(length(na.omit(MO2.NMS))))
)
OffPat.Group.All <- c()
for (i in 1:nrow(Challenge.Sum.All)) {
if (Challenge.Sum.All$Offspring.Acc[i] == 1 &
Challenge.Sum.All$Father.Acc[i] == 1) {
OffPat.Group.All[i] <- ("1")
} else if (
Challenge.Sum.All$Offspring.Acc[i] == 1 &
Challenge.Sum.All$Father.Acc[i] == 2) {
OffPat.Group.All[i] <- ("2")
} else if (
Challenge.Sum.All$Offspring.Acc[i] == 2 &
Challenge.Sum.All$Father.Acc[i] == 1) {
OffPat.Group.All[i] <- ("3")
} else if (
Challenge.Sum.All$Offspring.Acc[i] == 2 &
Challenge.Sum.All$Father.Acc[i] == 2) {
OffPat.Group.All[i] <- ("4")
}
}
Challenge.Sum.All$OffPat <- OffPat.Group.All
OffMat.Group.All <- c()
for (i in 1:nrow(Challenge.Sum.All)) {
if (Challenge.Sum.All$Offspring.Acc[i] == 1 &
Challenge.Sum.All$Mother.Acc[i] == 1) {
OffMat.Group.All[i] <- ("1")
} else if (
Challenge.Sum.All$Offspring.Acc[i] == 1 &
Challenge.Sum.All$Mother.Acc[i] == 2) {
OffMat.Group.All[i] <- ("2")
} else if (
Challenge.Sum.All$Offspring.Acc[i] == 2 &
Challenge.Sum.All$Mother.Acc[i] == 1) {
OffMat.Group.All[i] <- ("3")
} else if (
Challenge.Sum.All$Offspring.Acc[i] == 2 &
Challenge.Sum.All$Mother.Acc[i] == 2) {
OffMat.Group.All[i] <- ("4")
}
}
Challenge.Sum.All$OffMat <- OffMat.Group.All
MatPat.Group.All <- c()
for (i in 1:nrow(Challenge.Sum.All)) {
if (Challenge.Sum.All$Mother.Acc[i] == 1 &
Challenge.Sum.All$Father.Acc[i] == 1) {
MatPat.Group.All[i] <- ("1")
} else if (
Challenge.Sum.All$Mother.Acc[i] == 1 &
Challenge.Sum.All$Father.Acc[i] == 2) {
MatPat.Group.All[i] <- ("2")
} else if (
Challenge.Sum.All$Mother.Acc[i] == 2 &
Challenge.Sum.All$Father.Acc[i] == 1) {
MatPat.Group.All[i] <- ("3")
} else if (
Challenge.Sum.All$Mother.Acc[i] == 2 &
Challenge.Sum.All$Father.Acc[i] == 2) {
MatPat.Group.All[i] <- ("4")
}
}
Challenge.Sum.All$MatPat <- MatPat.Group.All
OffParError.DFacet <- ggplot(as.data.frame(Challenge.Sum.All), aes(
x = Challenge.Temp,
y = Mean, group = OffPat, shape = Father.Acc, colour = Offspring.Acc
)) +
theme_bw() +
my.theme +
geom_errorbar(
mapping = aes(
x = Challenge.Temp,
ymin = c(Lower.CI), ymax = c(Upper.CI)
), width = 0.2, size = 1,
position = position_dodge(width = 0.5)
) +
geom_point(size = 1.5, colour = "black", position = position_dodge(width = 0.5)) +
xlab("Challenge Temperature \n (°C)") +
ylab("Metabolic Rate (mg O2/h)") +
theme(legend.position = "top") +
scale_colour_manual(values = c("dodgerblue3", "maroon1")) +
guides(colour = guide_legend(title = "Offspring Acclimation \n Temperature")) +
guides(shape = guide_legend(title = "Paternal Acclimation \n Temperature")) +
theme(legend.title = element_text(size = 10))
## Facetting by maternal acclimation temperature first
OffParError.DFacet + facet_grid(rows = vars(Mother.Acc))
OffParError.SFacet <- ggplot(as.data.frame(Challenge.Sum.All), aes(
x = Challenge.Temp,
y = Mean, group = OffMat, shape = Mother.Acc, colour = Offspring.Acc
)) +
theme_bw() +
my.theme +
geom_errorbar(
mapping = aes(
x = Challenge.Temp,
ymin = c(Lower.CI), ymax = c(Upper.CI)
), width = 0.2, size = 1,
position = position_dodge(width = 0.5)
) +
geom_point(size = 1.5, colour = "black", position = position_dodge(width = 0.5)) +
xlab("Challenge Temperature \n (°C)") +
ylab("Metabolic Rate (mg O2/h)") +
theme(legend.position = "top") +
scale_colour_manual(values = c("slateblue", "gold2")) +
guides(colour = guide_legend(title = "Offspring Acclimation \n Temperature")) +
guides(shape = guide_legend(title = "Maternal Acclimation \n Temperature")) +
theme(legend.title = element_text(size = 10))
## Now facetting by paternal acclimation temperature
OffParError.SFacet + facet_grid(rows = vars(Father.Acc))
ParError.OFacet <- ggplot(as.data.frame(Challenge.Sum.All), aes(
x = Challenge.Temp,
y = Mean, group = MatPat, shape = Mother.Acc, colour = Father.Acc
)) +
theme_bw() +
my.theme +
geom_errorbar(
mapping = aes(
x = Challenge.Temp,
ymin = c(Lower.CI), ymax = c(Upper.CI)
), width = 0.2, size = 1,
position = position_dodge(width = 0.5)
) +
geom_point(size = 1.5, colour = "black", position = position_dodge(width = 0.5)) +
xlab("Challenge Temperature \n (°C)") +
ylab("Metabolic Rate (mg O2/h)") +
theme(legend.position = "top") +
scale_colour_manual(values = c("mediumseagreen", "mediumorchid")) +
guides(colour = guide_legend(title = "Paternal Acclimation \n Temperature")) +
guides(shape = guide_legend(title = "Maternal Acclimation \n Temperature")) +
theme(legend.title = element_text(size = 10))
## And lastly, facetting by offspring acclimation temperature
ParError.OFacet + facet_grid(rows = vars(Offspring.Acc))
Model construction
## Building cubic regression spline with low edf (here, 4) due to low curvature. Fitting spline to global data.
plot(MO2.NMS ~ Challenge.Temp, data = na.omit(BTChallenge))
curve.model <- gam(MO2.NMS ~ s(Challenge.Temp, bs = "cs", k = 7), na.action = na.exclude, data = BTChallenge)
## Results of gam
summary(curve.model)
##
## Family: gaussian
## Link function: identity
##
## Formula:
## MO2.NMS ~ s(Challenge.Temp, bs = "cs", k = 7)
##
## Parametric coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.776467 0.005796 134 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Approximate significance of smooth terms:
## edf Ref.df F p-value
## s(Challenge.Temp) 5.66 6 116.7 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## R-sq.(adj) = 0.213 Deviance explained = 21.5%
## GCV = 0.08691 Scale est. = 0.086686 n = 2580
plot(curve.model)
# Appears nicely fitting.
## Pulling out predicted values
smoother <- predict(curve.model, type = "response")
With.Fit <- BTChallenge
With.Fit$MO2.NMS <- smoother
With.Fit <- With.Fit %>%
mutate("Type" = "Fit")
Resp.Alone <- BTChallenge
Resp.Alone$Type <- "Response"
Total.Data <- rbind(With.Fit, Resp.Alone)
## Plotting the cubic regression spline
Total.Data %>%
ggplot(aes(x = Challenge.Temp, y = MO2.NMS, fill = Type, linetype = Type)) +
geom_point(pch = 21) +
scale_fill_manual(values = c("black", "magenta")) +
geom_smooth(method = "gam", formula = y ~ s(x, bs = "cr", k = 7))
## Again, appears fitting. Appending spline to dataframe
BTChallenge$bSpline <- smoother
# Constructing linear model
## Crossed random intercepts per parental ID, with crossed Offspring ID.
BTChallenge <- BTChallenge %>%
dplyr::rename("Degree.Day" = Degree.Dat)
BTChallenge$Degree.Day <- factor(BTChallenge$Degree.Day)
Dummy <- factor(1)
Dat <- groupedData(MO2.NMS ~ 1 | Dummy, BTChallenge)
Model1 <- lme(MO2.NMS ~ Mass + bSpline * Offspring.Acc * Father.Acc * Mother.Acc, random = pdBlocked(list(pdIdent(~ 0 + Father.ID), pdIdent(~ 0 + Mother.ID), pdIdent(~ 0 + Subject.ID), pdIdent(~ 0 + Degree.Day))), na.action = na.exclude, data = Dat)
# Assuming large degrees of freedom by estimating Z values (using glmmTMB):
require("glmmTMB")
Model1TMB <- glmmTMB(MO2.NMS ~ Mass + bSpline * Offspring.Acc * Father.Acc * Mother.Acc + (1 | Father.ID) + (1 | Mother.ID) + (1 | Subject.ID) + (1 | Degree.Day), na.action = na.exclude, data = BTChallenge)
# Diagnosing residuals
{
M1.Res.Alone <- ggplot(
data = data.frame("Res" = residuals(Model1, type = "normalized")),
aes(x = 1:nrow(Dat), y = Res)
) +
my.theme +
theme_bw() +
geom_point(size = 3, colour = "mediumseagreen", alpha = 0.5) +
xlab("Row Number") +
ylab("Deviance Residuals")
M1.Res.by.Fit <- ggplot(data = data.frame(
"Res" = residuals(Model1, type = "normalized"),
"Fit" = predict(Model1)
), aes(x = Fit, y = Res)) +
my.theme +
theme_bw() +
geom_point(size = 3, colour = "cornflowerblue", alpha = 0.5) +
xlab("Fitted Values") +
ylab("Deviance Residuals")
M1.Res.by.Resp <- ggplot(data = data.frame(
"Res" = residuals(Model1, type = "normalized"),
"Resp" = Dat$MO2.NMS
), aes(x = Resp, y = Res)) +
my.theme +
theme_bw() +
geom_point(size = 3, colour = "mediumorchid", alpha = 0.5) +
xlab("Metabolic Rate (mg O2/hour") +
ylab("Deviance Residuals")
M1.QQNorm.1 <- ggplot(
data.frame("Res" = residuals(Model1), "Fit" = predict(Model1)),
aes(sample = Res)
) +
stat_qq(colour = "gold") +
stat_qq_line() +
my.theme +
theme_bw()
}
grid.arrange(M1.Res.Alone, M1.Res.by.Fit, M1.Res.by.Resp, M1.QQNorm.1, nrow = 2, top = "Metabolic Rate Residuals")
# Large df
{
M1.Res.Alone.z <- ggplot(
data = data.frame("Res" = residuals(Model1TMB, type = "pearson")),
aes(x = 1:nrow(BTChallenge), y = Res)
) +
my.theme +
theme_bw() +
geom_point(size = 3, colour = "mediumseagreen", alpha = 0.5) +
xlab("Row Number") +
ylab("Pearson Residuals")
M1.Res.by.Fit.z <- ggplot(data = data.frame(
"Res" = residuals(Model1TMB, type = "pearson"),
"Fit" = predict(Model1)
), aes(x = Fit, y = Res)) +
my.theme +
theme_bw() +
geom_point(size = 3, colour = "cornflowerblue", alpha = 0.5) +
xlab("Fitted Values") +
ylab("Pearson Residuals")
M1.Res.by.Resp.z <- ggplot(data = data.frame(
"Res" = residuals(Model1TMB, type = "pearson"),
"Resp" = BTChallenge$MO2.NMS
), aes(x = Resp, y = Res)) +
my.theme +
theme_bw() +
geom_point(size = 3, colour = "mediumorchid", alpha = 0.5) +
xlab("Metabolic Rate (mg O2/hour") +
ylab("Pearson Residuals")
M1.QQNorm.1.z <- ggplot(
data.frame("Res" = residuals(Model1TMB), "Fit" = predict(Model1TMB)),
aes(sample = Res)
) +
stat_qq(colour = "gold") +
stat_qq_line() +
my.theme +
theme_bw()
}
grid.arrange(M1.Res.Alone.z, M1.Res.by.Fit.z, M1.Res.by.Resp.z, M1.QQNorm.1.z, nrow = 2, top = "Metabolic Rate Residuals")
## T model: Straying from normality in upper quantiles. Fan notable in residuals by fitteds. Perhaps weight by challenge temperature. First, testing for autocorrelation.
## Z model: Bias in residuals. Very well could be a consequence of autocorrelation.
ggplot(with(
acf(na.omit(residuals(Model1, type = "normalized")), plot = "FALSE"),
data.frame(lag, acf)
), aes(x = lag, y = acf)) +
geom_bar(
stat = "identity",
position = "identity", fill = "slateblue", colour = "black"
) +
theme_bw() +
my.theme.diag +
ylab("Autocorrelation \n Factor")
rho <- acf(na.omit(residuals(Model1, type = "normalized")), plot = F)$acf[2]
print(rho)
## [1] 0.2207105
ggplot(with(
acf(na.omit(residuals(Model1TMB, type = "pearson")), plot = "FALSE"),
data.frame(lag, acf)
), aes(x = lag, y = acf)) +
geom_bar(
stat = "identity",
position = "identity", fill = "slateblue", colour = "black"
) +
theme_bw() +
my.theme.diag +
ylab("Autocorrelation \n Factor")
rho <- acf(na.omit(residuals(Model1, type = "normalized")), plot = F)$acf[2]
print(rho)
## [1] 0.2207105
# Minor autocorrelation between measures in both models model. Estimating rho.
# Faceting by paternal acclimation temperature for clearer viewing.
Dat %>%
mutate("Res" = residuals(Model1, type = "normalized")) %>%
ggplot(aes(x = 1:nrow(Dat), y = Res, fill = Subject.ID, colour = Subject.ID)) +
my.theme +
theme_bw() +
geom_point(
size = 1, colour = "mediumseagreen",
alpha = 0.5, pch = 21
) +
xlab("Row Number") +
ylab("Deviance Residuals") +
theme(legend.position = "NULL") +
geom_line(alpha = 0.4) +
facet_grid(
rows =
Dat$Father.Acc
)
Dat %>%
mutate("Res" = residuals(Model1TMB, type = "pearson")) %>%
ggplot(aes(x = 1:nrow(Dat), y = Res, fill = Subject.ID, colour = Subject.ID)) +
my.theme +
theme_bw() +
geom_point(
size = 1, colour = "mediumseagreen",
alpha = 0.5, pch = 21
) +
xlab("Row Number") +
ylab("Pearson Residuals") +
theme(legend.position = "NULL") +
geom_line(alpha = 0.4) +
facet_grid(
rows =
Dat$Father.Acc
)
## Similar lag to models in Penney et al (In Press). Residuals, however, appear evenly spread.
# Adjusting for autocorrelation.
Model.AR <- lme(MO2.NMS ~ Mass + bSpline * Offspring.Acc * Father.Acc * Mother.Acc, random = pdBlocked(list(pdIdent(~ 0 + Father.ID), pdIdent(~ 0 + Mother.ID), pdIdent(~ 0 + Subject.ID), pdIdent(~ 0 + Degree.Day))), correlation = corAR1(), na.action = na.exclude, data = Dat)
BTChallenge$AR_Group <- factor(with(BTChallenge, paste(Father.ID, Mother.ID, Subject.ID, sep = "_")))
BTChallenge <- BTChallenge %>%
mutate("Time" = Challenge.Temp - 14) %>%
mutate("Time" = factor(Time)) %>%
arrange(AR_Group, Time) %>%
as.data.frame(.)
ModelTMB.AR <- glmmTMB(MO2.NMS ~ Mass + bSpline * Offspring.Acc * Father.Acc * Mother.Acc + (1 | Father.ID) + (1 | Mother.ID) + (1 | Subject.ID) + (1 | Degree.Day) + ar1(Time + 0 | AR_Group), na.action = na.exclude, data = BTChallenge)
## Diagnostics again.
# T-model
{
MAR.Res.Alone <- ggplot(data = data.frame("Res" = residuals(Model.AR,
type =
"normalized"
)), aes(x = 1:nrow(Dat), y = Res)) +
my.theme +
theme_bw() +
geom_point(size = 3, colour = "mediumseagreen", alpha = 0.5) +
xlab("Row Number") +
ylab("Deviance Residuals")
MAR.Res.by.Fit <- ggplot(data = data.frame("Res" = residuals(Model.AR,
type =
"normalized"
), "Fit" = predict(Model.AR)), aes(x = Fit, y = Res)) +
my.theme +
theme_bw() +
geom_point(size = 3, colour = "cornflowerblue", alpha = 0.5) +
xlab("Fitted Values") +
ylab("Deviance Residuals")
MAR.Res.by.Resp <- ggplot(data = data.frame("Res" = residuals(Model.AR,
type =
"normalized"
), "Resp" = Dat$MO2.NMS), aes(x = Resp, y = Res)) +
my.theme +
theme_bw() +
geom_point(size = 3, colour = "mediumorchid", alpha = 0.5) +
xlab("Metabolic Rate (mg O2/hour") +
ylab("Deviance Residuals")
MAR.QQNorm.1 <- ggplot(
data.frame("Res" = residuals(Model.AR), "Fit" = predict(Model.AR)),
aes(sample = Res)
) +
stat_qq(colour = "gold") +
stat_qq_line() +
my.theme +
theme_bw()
}
grid.arrange(MAR.Res.Alone, MAR.Res.by.Fit, MAR.Res.by.Resp, MAR.QQNorm.1, nrow = 2, top = "Metabolic Rate Residuals - Autoregressive Correlation Structure")
# Z-model
{
MAR.Res.Alone.z <- ggplot(data = data.frame("Res" = residuals(ModelTMB.AR,
type =
"response"
)), aes(x = 1:nrow(BTChallenge), y = Res)) +
my.theme +
theme_bw() +
geom_point(size = 3, colour = "mediumseagreen", alpha = 0.5) +
xlab("Row Number") +
ylab("Response Residuals")
MAR.Res.by.Fit.z <- ggplot(data = data.frame("Res" = residuals(ModelTMB.AR,
type =
"response"
), "Fit" = predict(ModelTMB.AR)), aes(x = Fit, y = Res)) +
my.theme +
theme_bw() +
geom_point(size = 3, colour = "cornflowerblue", alpha = 0.5) +
xlab("Fitted Values") +
ylab("Response Residuals")
MAR.Res.by.Resp.z <- ggplot(data = data.frame("Res" = residuals(ModelTMB.AR,
type =
"response"
), "Resp" = BTChallenge$MO2.NMS), aes(x = Resp, y = Res)) +
my.theme +
theme_bw() +
geom_point(size = 3, colour = "mediumorchid", alpha = 0.5) +
xlab("Metabolic Rate (mg O2/hour") +
ylab("Response Residuals")
MAR.QQNorm.1.z <- ggplot(
data.frame("Res" = residuals(ModelTMB.AR), "Fit" = predict(ModelTMB.AR)),
aes(sample = Res)
) +
stat_qq(colour = "gold") +
stat_qq_line() +
my.theme +
theme_bw()
}
grid.arrange(MAR.Res.Alone.z, MAR.Res.by.Fit.z, MAR.Res.by.Resp.z, MAR.QQNorm.1.z, nrow = 2, top = "Metabolic Rate Residuals - Autoregressive Correlation Structure")
# Very similar outcome to previous models. In glmmTMB model, similarities could merely be a consequence of plotted residuals not correcting for autocorrelation. Checking
plot(residuals(ModelTMB.AR, type = "response") ~ residuals(Model1TMB, type = "response"))
plot(residuals(ModelTMB.AR, type = "pearson") ~ residuals(Model1TMB, type = "pearson"))
# Yes, this appears to be the case. Using DHARMa to assess residuals
require("DHARMa")
BTChallenge_Clean <- na.omit(BTChallenge)
BTChallenge_Clean$Time_Int <- as.integer(BTChallenge_Clean$Time)
ModelTMB.AR <- glmmTMB(MO2.NMS ~ Mass + bSpline * Offspring.Acc * Father.Acc * Mother.Acc + (1 | Father.ID) + (1 | Mother.ID) + (1 | Subject.ID) + (1 | Degree.Day) + ar1(Time + 0 | AR_Group), na.action = na.omit, data = BTChallenge_Clean)
dharma_resid <- simulateResiduals(fittedModel = ModelTMB.AR)
plotResiduals(dharma_resid, quantreg = T) # Yes, residuals already appear cleaner.
dharma_residAR <- recalculateResiduals(dharma_resid, group = BTChallenge_Clean$Time)
testTemporalAutocorrelation(simulationOutput = dharma_residAR, time = sort(unique(BTChallenge_Clean$Time)))
##
## Durbin-Watson test
##
## data: simulationOutput$scaledResiduals ~ 1
## DW = 1.7116, p-value = 0.5655
## alternative hypothesis: true autocorrelation is not 0
ggplot(with(
acf(na.omit(residuals(Model.AR, type = "normalized")), plot = "FALSE"),
data.frame(lag, acf)
), aes(x = lag, y = acf)) +
geom_bar(
stat = "identity",
position = "identity", fill = "slateblue", colour = "black"
) +
theme_bw() +
my.theme.diag +
ylab("Autocorrelation \n Factor")
## Yes, autocorrelation and residuals do appear okay in both models.
# Again, faceting by paternal acclimation temperature for clearer viewing.
Dat %>%
mutate("Res" = residuals(Model.AR, type = "normalized")) %>%
ggplot(aes(x = 1:nrow(Dat), y = Res, fill = Subject.ID, colour = Subject.ID)) +
my.theme +
theme_bw() +
geom_point(
size = 1, colour = "mediumseagreen",
alpha = 0.5, pch = 21
) +
xlab("Row Number") +
ylab("Deviance Residuals") +
theme(legend.position = "NULL") +
geom_line(alpha = 0.4) +
facet_grid(
rows =
Dat$Father.Acc
)
BTChallenge_Clean %>%
mutate("Res" = dharma_resid$scaledResiduals) %>%
ggplot(aes(x = 1:nrow(BTChallenge_Clean), y = Res, fill = Subject.ID, colour = Subject.ID)) +
my.theme +
theme_bw() +
geom_point(
size = 1, colour = "mediumseagreen",
alpha = 0.5, pch = 21
) +
xlab("Row Number") +
ylab("Deviance Residuals") +
theme(legend.position = "NULL") +
geom_line(alpha = 0.4) +
facet_grid(
rows =
BTChallenge_Clean$Father.Acc
)
# Good.
## Addressing tail in residuals
Dat %>%
mutate("Res" = residuals(Model1, type = "normalized"), "Fit" = predict(Model.AR)) %>%
ggplot(aes(sample = Res, fill = Challenge.Temp)) +
geom_qq(pch = 21, size = 3, alpha = 0.5, aes(fill = factor(Dat$Challenge.Temp))) +
stat_qq_line() +
my.theme +
theme_bw()
G_Col <- c(Cold = "dodgerblue3", Warm = "gold2", Hot = "firebrick")
Dat %>%
mutate("Res" = abs(residuals(Model1, type = "normalized"))) %>%
mutate("Temp_Group" = ifelse(Challenge.Temp < 20, "Cold",
ifelse(Challenge.Temp > 25, "Hot", "Warm")
)) %>%
mutate("Temp_Group" = factor(Temp_Group, levels = c("Cold", "Warm", "Hot"))) %>%
ggplot(aes(x = Res, fill = Temp_Group)) +
geom_density(alpha = 0.4, colour = "black", adjust = 2) +
scale_fill_manual(values = G_Col) +
theme_classic()
BTChallenge_Clean %>%
mutate("Res" = dharma_resid$scaledResiduals) %>%
mutate("Temp_Group" = ifelse(Challenge.Temp < 20, "Cold",
ifelse(Challenge.Temp > 25, "Hot", "Warm")
)) %>%
mutate("Temp_Group" = factor(Temp_Group, levels = c("Cold", "Warm", "Hot"))) %>%
ggplot(aes(x = Res, fill = Temp_Group)) +
geom_density(alpha = 0.4, colour = "black", adjust = 2) +
scale_fill_manual(values = G_Col) +
theme_classic()
# Slightly disproportionate pulling in cold and warm groups. Weighting by challenge temperature
Model.AR2 <- lme(MO2.NMS ~ Mass + bSpline * Offspring.Acc * Father.Acc * Mother.Acc, random = pdBlocked(list(pdIdent(~ 0 + Father.ID), pdIdent(~ 0 + Mother.ID), pdIdent(~ 0 + Subject.ID), pdIdent(~ 0 + Degree.Day))), correlation = corAR1(), weights = ~ 1 / Challenge.Temp, na.action = na.exclude, data = Dat)
ModelTMB2.AR <- glmmTMB(MO2.NMS ~ Mass + bSpline * Offspring.Acc * Father.Acc * Mother.Acc + (1 | Father.ID) + (1 | Mother.ID) + (1 | Subject.ID) + (1 | Degree.Day) + ar1(Time + 0 | AR_Group), weights = 1 / Challenge.Temp, na.action = na.omit, data = BTChallenge_Clean)
# glmmTMB model failing to converge. Assessing residuals of lme model.
{
MAR2.Res.Alone <- ggplot(data = data.frame("Res" = residuals(Model.AR2,
type =
"normalized"
)), aes(x = 1:nrow(Dat), y = Res)) +
my.theme +
theme_bw() +
geom_point(size = 3, colour = "mediumseagreen", alpha = 0.5) +
xlab("Row Number") +
ylab("Deviance Residuals")
MAR2.Res.by.Fit <- ggplot(data = data.frame("Res" = residuals(Model.AR2,
type =
"normalized"
), "Fit" = predict(Model.AR2)), aes(x = Fit, y = Res)) +
my.theme +
theme_bw() +
geom_point(size = 3, colour = "cornflowerblue", alpha = 0.5) +
xlab("Fitted Values") +
ylab("Deviance Residuals")
MAR2.Res.by.Resp <- ggplot(data = data.frame("Res" = residuals(Model.AR2,
type =
"normalized"
), "Resp" = Dat$MO2.NMS), aes(x = Resp, y = Res)) +
my.theme +
theme_bw() +
geom_point(size = 3, colour = "mediumorchid", alpha = 0.5) +
xlab("Metabolic Rate (mg O2/hour") +
ylab("Deviance Residuals")
MAR2.QQNorm <- ggplot(data.frame("Res" = residuals(Model.AR2, type = "normalized"), "Fit" = predict(Model.AR2)), aes(sample = Res)) +
stat_qq(colour = "gold") +
stat_qq_line() +
my.theme +
theme_bw()
}
grid.arrange(MAR2.Res.Alone, MAR2.Res.by.Fit, MAR2.Res.by.Resp, MAR2.QQNorm, nrow = 2, top = "Metabolic Rate Residuals - Autoregressive Correlation Structure and Weighted by Challenge Temp.")
## Tail persisting, but subtly. Plotting histogram of residuals.
Dat %>%
mutate("Res" = residuals(Model.AR2, type = "normalized")) %>%
ggplot(aes(x = Res)) +
geom_histogram(fill = "magenta", colour = "black", size = 1, bins = 20) +
xlab("Normalized Residuals")
## Appears quite nice, though a few notable outliers. Identifying outliers.
Dat <- Dat %>%
mutate("Res" = residuals(Model.AR2, type = "normalized"))
Grouped <- rbind(subset(Dat, Res < with(na.omit(Dat), mean(Res) - 3 * sd(Res))) %>%
mutate("Res.Group" = "Low.Out"), subset(Dat, Res > with(na.omit(Dat), mean(Res) + 3 * sd(Res))) %>%
mutate("Res.Group" = "High.Out"), subset(Dat, Res >= with(na.omit(Dat), mean(Res) - 3 * sd(Res)) & Res <= with(na.omit(Dat), mean(Res) + 3 * sd(Res))) %>%
mutate("Res.Group" = "Normal"))
G1 <- data.frame(subset(Grouped, Res.Group == "High.Out") %>%
count(Offspring.Acc), "Group" = "Offspring.Acc")
colnames(G1) <- c("Group.Number", "Count", "Group")
G2 <- data.frame(subset(Grouped, Res.Group == "High.Out") %>%
count(Father.Acc), "Group" = "Father.Acc")
colnames(G2) <- c("Group.Number", "Count", "Group")
G3 <- data.frame(subset(Grouped, Res.Group == "High.Out") %>%
count(Mother.Acc), "Group" = "Mother.Acc")
colnames(G3) <- c("Group.Number", "Count", "Group")
rbind(G1, G2, G3)
## Group.Number Count Group
## 1 1 5 Offspring.Acc
## 2 2 10 Offspring.Acc
## 3 1 11 Father.Acc
## 4 2 4 Father.Acc
## 5 1 5 Mother.Acc
## 6 2 10 Mother.Acc
## Possible interactive assortment. Checking below.
subset(Grouped, Res.Group == "High.Out") %>%
count(Father.Acc:Mother.Acc)
## Grouped Data: MO2.NMS ~ 1 | Dummy
## Father.Acc:Mother.Acc n
## 1 1:1 3
## 2 1:2 8
## 3 2:1 2
## 4 2:2 2
## Not distinct. Checking trends accross acclimation temperature
subset(Grouped, Res.Group == "High.Out") %>%
count(Challenge.Temp)
## Grouped Data: MO2.NMS ~ 1 | Dummy
## Challenge.Temp n
## 1 18 1
## 2 19 2
## 3 20 1
## 4 21 1
## 5 23 1
## 6 25 1
## 7 26 3
## 8 27 3
## 9 28 2
## Not clear either. Residual trends do not appear systematic. Continuing with previous model as final.
Plotting Predicted Values
Dat$Pred <- predict(Model.AR)
Challenge.Sum.Pred <- Dat %>%
group_by(Challenge.Temp, Offspring.Acc, Father.Acc, Mother.Acc) %>%
summarize(
Mean = mean(Pred, na.rm = TRUE), Upper.CI = mean(Pred,
na.rm = TRUE
) + 1.96 * (sd(Pred, na.rm = TRUE) / sqrt(length(na.omit(Pred)))),
Lower.CI = mean(Pred, na.rm = TRUE) - 1.96 * (sd(Pred, na.rm = TRUE) /
sqrt(length(na.omit(Pred))))
)
OffPat.Group.Pred <- c()
for (i in 1:nrow(Challenge.Sum.Pred)) {
if (Challenge.Sum.Pred$Offspring.Acc[i] == 1 &
Challenge.Sum.Pred$Father.Acc[i] == 1) {
OffPat.Group.Pred[i] <- ("1")
} else if (
Challenge.Sum.Pred$Offspring.Acc[i] == 1 &
Challenge.Sum.Pred$Father.Acc[i] == 2) {
OffPat.Group.Pred[i] <- ("2")
} else if (
Challenge.Sum.Pred$Offspring.Acc[i] == 2 &
Challenge.Sum.Pred$Father.Acc[i] == 1) {
OffPat.Group.Pred[i] <- ("3")
} else if (
Challenge.Sum.Pred$Offspring.Acc[i] == 2 &
Challenge.Sum.Pred$Father.Acc[i] == 2) {
OffPat.Group.Pred[i] <- ("4")
}
}
Challenge.Sum.Pred$OffPat <- OffPat.Group.Pred
# Faceting by maternal acclimation temperature.
PredError <- ggplot(as.data.frame(Challenge.Sum.Pred), aes(
x = Challenge.Temp,
y = Mean, group = OffPat, shape = Father.Acc, colour = Offspring.Acc
)) +
theme_bw() +
my.theme +
geom_errorbar(
mapping = aes(
x = Challenge.Temp,
ymin = c(Lower.CI), ymax = c(Upper.CI)
), width = 0.2, size = 1,
position = position_dodge(width = 0.5)
) +
geom_point(size = 1.5, colour = "black", position = position_dodge(width = 0.5)) +
xlab("Challenge Temperature \n (°C)") +
ylab("Metabolic Rate (mg O2/h)") +
theme(legend.position = "top") +
scale_colour_manual(values = c("dodgerblue3", "maroon1")) +
guides(colour = guide_legend(title = "Offspring Acclimation \n Temperature")) +
guides(shape = guide_legend(title = "Paternal Acclimation \n Temperature")) +
theme(legend.title = element_text(size = 10))
PredError + facet_grid(rows = vars(Mother.Acc))
OffMat.Group.Pred <- c()
for (i in 1:nrow(Challenge.Sum.Pred)) {
if (Challenge.Sum.Pred$Offspring.Acc[i] == 1 &
Challenge.Sum.Pred$Mother.Acc[i] == 1) {
OffMat.Group.Pred[i] <- ("1")
} else if (
Challenge.Sum.Pred$Offspring.Acc[i] == 1 &
Challenge.Sum.Pred$Mother.Acc[i] == 2) {
OffMat.Group.Pred[i] <- ("2")
} else if (
Challenge.Sum.Pred$Offspring.Acc[i] == 2 &
Challenge.Sum.Pred$Mother.Acc[i] == 1) {
OffMat.Group.Pred[i] <- ("3")
} else if (
Challenge.Sum.Pred$Offspring.Acc[i] == 2 &
Challenge.Sum.Pred$Mother.Acc[i] == 2) {
OffMat.Group.Pred[i] <- ("4")
}
}
Challenge.Sum.Pred$OffMat <- OffMat.Group.Pred
# And faceting by paternal acclimation temperature
PredError <- ggplot(as.data.frame(Challenge.Sum.Pred), aes(
x = Challenge.Temp,
y = Mean, group = OffMat, shape = Mother.Acc, colour = Offspring.Acc
)) +
theme_bw() +
my.theme +
geom_errorbar(
mapping = aes(
x = Challenge.Temp,
ymin = c(Lower.CI), ymax = c(Upper.CI)
), width = 0.2, size = 1,
position = position_dodge(width = 0.5)
) +
geom_point(size = 1.5, colour = "black", position = position_dodge(width = 0.5)) +
xlab("Challenge Temperature \n (°C)") +
ylab("Metabolic Rate (mg O2/h)") +
theme(legend.position = "top") +
scale_colour_manual(values = c("cornflowerblue", "gold2")) +
guides(colour = guide_legend(title = "Offspring Acclimation \n Temperature")) +
guides(shape = guide_legend(title = "Maternal Acclimation \n Temperature")) +
theme(legend.title = element_text(size = 10))
PredError + facet_grid(rows = vars(Father.Acc))
## Striking model fit. Checking glmmTMB model
BTChallenge_Clean$Pred <- predict(ModelTMB.AR)
Challenge.Sum.Pred <- BTChallenge_Clean %>%
group_by(Challenge.Temp, Offspring.Acc, Father.Acc, Mother.Acc) %>%
summarize(
Mean = mean(Pred, na.rm = TRUE), Upper.CI = mean(Pred,
na.rm = TRUE
) + 1.96 * (sd(Pred, na.rm = TRUE) / sqrt(length(na.omit(Pred)))),
Lower.CI = mean(Pred, na.rm = TRUE) - 1.96 * (sd(Pred, na.rm = TRUE) /
sqrt(length(na.omit(Pred))))
)
OffPat.Group.Pred <- c()
for (i in 1:nrow(Challenge.Sum.Pred)) {
if (Challenge.Sum.Pred$Offspring.Acc[i] == 1 &
Challenge.Sum.Pred$Father.Acc[i] == 1) {
OffPat.Group.Pred[i] <- ("1")
} else if (
Challenge.Sum.Pred$Offspring.Acc[i] == 1 &
Challenge.Sum.Pred$Father.Acc[i] == 2) {
OffPat.Group.Pred[i] <- ("2")
} else if (
Challenge.Sum.Pred$Offspring.Acc[i] == 2 &
Challenge.Sum.Pred$Father.Acc[i] == 1) {
OffPat.Group.Pred[i] <- ("3")
} else if (
Challenge.Sum.Pred$Offspring.Acc[i] == 2 &
Challenge.Sum.Pred$Father.Acc[i] == 2) {
OffPat.Group.Pred[i] <- ("4")
}
}
Challenge.Sum.Pred$OffPat <- OffPat.Group.Pred
# Faceting by maternal acclimation temperature again.
PredError <- ggplot(as.data.frame(Challenge.Sum.Pred), aes(
x = Challenge.Temp,
y = Mean, group = OffPat, shape = Father.Acc, colour = Offspring.Acc
)) +
theme_bw() +
my.theme +
geom_errorbar(
mapping = aes(
x = Challenge.Temp,
ymin = c(Lower.CI), ymax = c(Upper.CI)
), width = 0.2, size = 1,
position = position_dodge(width = 0.5)
) +
geom_point(size = 1.5, colour = "black", position = position_dodge(width = 0.5)) +
xlab("Challenge Temperature \n (°C)") +
ylab("Metabolic Rate (mg O2/h)") +
theme(legend.position = "top") +
scale_colour_manual(values = c("dodgerblue3", "maroon1")) +
guides(colour = guide_legend(title = "Offspring Acclimation \n Temperature")) +
guides(shape = guide_legend(title = "Paternal Acclimation \n Temperature")) +
theme(legend.title = element_text(size = 10))
PredError + facet_grid(rows = vars(Mother.Acc))
OffMat.Group.Pred <- c()
for (i in 1:nrow(Challenge.Sum.Pred)) {
if (Challenge.Sum.Pred$Offspring.Acc[i] == 1 &
Challenge.Sum.Pred$Mother.Acc[i] == 1) {
OffMat.Group.Pred[i] <- ("1")
} else if (
Challenge.Sum.Pred$Offspring.Acc[i] == 1 &
Challenge.Sum.Pred$Mother.Acc[i] == 2) {
OffMat.Group.Pred[i] <- ("2")
} else if (
Challenge.Sum.Pred$Offspring.Acc[i] == 2 &
Challenge.Sum.Pred$Mother.Acc[i] == 1) {
OffMat.Group.Pred[i] <- ("3")
} else if (
Challenge.Sum.Pred$Offspring.Acc[i] == 2 &
Challenge.Sum.Pred$Mother.Acc[i] == 2) {
OffMat.Group.Pred[i] <- ("4")
}
}
Challenge.Sum.Pred$OffMat <- OffMat.Group.Pred
# Paternal acclimation temperature?
PredError <- ggplot(as.data.frame(Challenge.Sum.Pred), aes(
x = Challenge.Temp,
y = Mean, group = OffMat, shape = Mother.Acc, colour = Offspring.Acc
)) +
theme_bw() +
my.theme +
geom_errorbar(
mapping = aes(
x = Challenge.Temp,
ymin = c(Lower.CI), ymax = c(Upper.CI)
), width = 0.2, size = 1,
position = position_dodge(width = 0.5)
) +
geom_point(size = 1.5, colour = "black", position = position_dodge(width = 0.5)) +
xlab("Challenge Temperature \n (°C)") +
ylab("Metabolic Rate (mg O2/h)") +
theme(legend.position = "top") +
scale_colour_manual(values = c("cornflowerblue", "gold2")) +
guides(colour = guide_legend(title = "Offspring Acclimation \n Temperature")) +
guides(shape = guide_legend(title = "Maternal Acclimation \n Temperature")) +
theme(legend.title = element_text(size = 10))
PredError + facet_grid(rows = vars(Father.Acc))
# Pulling out degrees of freedom using Satterthwaite method (in lmerTest)
DF.Model <- lmer(MO2.NMS ~ Mass + bSpline * Offspring.Acc * Father.Acc * Mother.Acc + (1 | Father.ID) + (1 | Mother.ID) + (1 | Subject.ID) + (1 | Degree.Day), na.action = na.exclude, data = Dat)
DF <- data.frame(summary(DF.Model)$coefficients)$df
t_vals <- data.frame(summary(Model.AR)$tTab)$t.value
p_vals <- c()
for (i in 1:length(t_vals)) {
if (t_vals[i] > 0) {
p_vals[i] <- 2 * (pt(t_vals[i], df = DF[i], lower = FALSE))
} else if (t_vals[i] < 0) {
p_vals[i] <- 2 * (pt(t_vals[i], df = DF[i], lower = TRUE))
}
}
coef_output <- data.frame(summary(Model.AR)$tTab)
coef_output$p.value <- p_vals
coef_output$DF <- DF
# Printing full model output
summary(Model.AR)
## Linear mixed-effects model fit by REML
## Data: Dat
## AIC BIC logLik
## -797.621 -663.1764 421.8105
##
## Random effects:
## Composite Structure: Blocked
##
## Block 1: Father.ID1, Father.ID2, Father.ID3, Father.ID4
## Formula: ~0 + Father.ID | Dummy
## Structure: Multiple of an Identity
## Father.ID1 Father.ID2 Father.ID3 Father.ID4
## StdDev: 0.07769035 0.07769035 0.07769035 0.07769035
##
## Block 2: Mother.ID1, Mother.ID2, Mother.ID3, Mother.ID4
## Formula: ~0 + Mother.ID | Dummy
## Structure: Multiple of an Identity
## Mother.ID1 Mother.ID2 Mother.ID3 Mother.ID4
## StdDev: 0.005968334 0.005968334 0.005968334 0.005968334
##
## Block 3: Subject.ID1, Subject.ID2, Subject.ID3, Subject.ID4, Subject.ID5, Subject.ID6, Subject.ID7, Subject.ID8, Subject.ID9, Subject.ID10, Subject.ID11, Subject.ID12, Subject.ID13, Subject.ID14, Subject.ID15, Subject.ID16, Subject.ID17, Subject.ID18, Subject.ID19, Subject.ID20, Subject.ID21, Subject.ID22, Subject.ID23, Subject.ID24, Subject.ID25, Subject.ID26, Subject.ID27, Subject.ID28, Subject.ID29, Subject.ID30, Subject.ID31, Subject.ID32, Subject.ID33, Subject.ID34, Subject.ID35, Subject.ID36, Subject.ID37, Subject.ID38, Subject.ID39, Subject.ID40, Subject.ID41, Subject.ID42, Subject.ID43, Subject.ID44, Subject.ID45, Subject.ID46, Subject.ID47, Subject.ID48, Subject.ID49, Subject.ID50, Subject.ID51, Subject.ID52, Subject.ID53, Subject.ID54, Subject.ID55, Subject.ID56, Subject.ID57, Subject.ID58, Subject.ID59, Subject.ID60, Subject.ID61, Subject.ID62, Subject.ID63, Subject.ID64, Subject.ID65, Subject.ID66, Subject.ID67, Subject.ID68, Subject.ID69, Subject.ID70, Subject.ID71, Subject.ID72, Subject.ID73, Subject.ID74, Subject.ID75, Subject.ID76, Subject.ID77, Subject.ID78, Subject.ID79, Subject.ID80, Subject.ID81, Subject.ID82, Subject.ID83, Subject.ID84, Subject.ID85, Subject.ID86, Subject.ID87, Subject.ID88, Subject.ID89, Subject.ID90, Subject.ID91, Subject.ID92, Subject.ID93, Subject.ID94, Subject.ID95, Subject.ID96, Subject.ID97, Subject.ID98, Subject.ID99, Subject.ID100, Subject.ID101, Subject.ID102, Subject.ID103, Subject.ID104, Subject.ID105, Subject.ID106, Subject.ID107, Subject.ID108, Subject.ID109, Subject.ID110, Subject.ID111, Subject.ID112, Subject.ID113, Subject.ID114, Subject.ID115, Subject.ID116, Subject.ID117, Subject.ID118, Subject.ID119, Subject.ID120, Subject.ID121, Subject.ID122, Subject.ID123, Subject.ID124, Subject.ID125, Subject.ID126, Subject.ID127, Subject.ID128, Subject.ID129, Subject.ID130, Subject.ID131, Subject.ID132, Subject.ID133, Subject.ID134, Subject.ID135, Subject.ID136, Subject.ID137, Subject.ID138, Subject.ID139, Subject.ID140, Subject.ID141, Subject.ID142, Subject.ID143, Subject.ID144, Subject.ID145, Subject.ID146, Subject.ID147, Subject.ID148, Subject.ID149, Subject.ID150, Subject.ID151, Subject.ID152, Subject.ID153, Subject.ID154, Subject.ID155, Subject.ID156, Subject.ID157, Subject.ID158, Subject.ID159, Subject.ID160, Subject.ID161, Subject.ID162, Subject.ID163, Subject.ID164, Subject.ID165, Subject.ID166, Subject.ID167, Subject.ID168, Subject.ID169, Subject.ID170, Subject.ID171, Subject.ID172, Subject.ID173, Subject.ID174, Subject.ID175, Subject.ID176, Subject.ID177, Subject.ID178, Subject.ID179, Subject.ID180, Subject.ID181, Subject.ID182, Subject.ID183, Subject.ID184, Subject.ID185, Subject.ID186, Subject.ID187, Subject.ID188, Subject.ID189, Subject.ID190, Subject.ID191, Subject.ID192, Subject.ID193, Subject.ID194, Subject.ID195, Subject.ID196, Subject.ID197, Subject.ID198, Subject.ID199, Subject.ID200, Subject.ID201, Subject.ID202, Subject.ID203, Subject.ID204, Subject.ID205, Subject.ID206, Subject.ID207, Subject.ID208, Subject.ID209, Subject.ID210, Subject.ID211, Subject.ID212, Subject.ID213, Subject.ID214, Subject.ID215, Subject.ID216, Subject.ID217, Subject.ID218, Subject.ID219, Subject.ID220, Subject.ID221, Subject.ID222, Subject.ID223, Subject.ID224, Subject.ID225, Subject.ID227, Subject.ID228, Subject.ID229, Subject.ID230
## Formula: ~0 + Subject.ID | Dummy
## Structure: Multiple of an Identity
## Subject.ID1 Subject.ID2 Subject.ID3 Subject.ID4 Subject.ID5 Subject.ID6
## StdDev: 0.1471092 0.1471092 0.1471092 0.1471092 0.1471092 0.1471092
## Subject.ID7 Subject.ID8 Subject.ID9 Subject.ID10 Subject.ID11
## StdDev: 0.1471092 0.1471092 0.1471092 0.1471092 0.1471092
## Subject.ID12 Subject.ID13 Subject.ID14 Subject.ID15 Subject.ID16
## StdDev: 0.1471092 0.1471092 0.1471092 0.1471092 0.1471092
## Subject.ID17 Subject.ID18 Subject.ID19 Subject.ID20 Subject.ID21
## StdDev: 0.1471092 0.1471092 0.1471092 0.1471092 0.1471092
## Subject.ID22 Subject.ID23 Subject.ID24 Subject.ID25 Subject.ID26
## StdDev: 0.1471092 0.1471092 0.1471092 0.1471092 0.1471092
## Subject.ID27 Subject.ID28 Subject.ID29 Subject.ID30 Subject.ID31
## StdDev: 0.1471092 0.1471092 0.1471092 0.1471092 0.1471092
## Subject.ID32 Subject.ID33 Subject.ID34 Subject.ID35 Subject.ID36
## StdDev: 0.1471092 0.1471092 0.1471092 0.1471092 0.1471092
## Subject.ID37 Subject.ID38 Subject.ID39 Subject.ID40 Subject.ID41
## StdDev: 0.1471092 0.1471092 0.1471092 0.1471092 0.1471092
## Subject.ID42 Subject.ID43 Subject.ID44 Subject.ID45 Subject.ID46
## StdDev: 0.1471092 0.1471092 0.1471092 0.1471092 0.1471092
## Subject.ID47 Subject.ID48 Subject.ID49 Subject.ID50 Subject.ID51
## StdDev: 0.1471092 0.1471092 0.1471092 0.1471092 0.1471092
## Subject.ID52 Subject.ID53 Subject.ID54 Subject.ID55 Subject.ID56
## StdDev: 0.1471092 0.1471092 0.1471092 0.1471092 0.1471092
## Subject.ID57 Subject.ID58 Subject.ID59 Subject.ID60 Subject.ID61
## StdDev: 0.1471092 0.1471092 0.1471092 0.1471092 0.1471092
## Subject.ID62 Subject.ID63 Subject.ID64 Subject.ID65 Subject.ID66
## StdDev: 0.1471092 0.1471092 0.1471092 0.1471092 0.1471092
## Subject.ID67 Subject.ID68 Subject.ID69 Subject.ID70 Subject.ID71
## StdDev: 0.1471092 0.1471092 0.1471092 0.1471092 0.1471092
## Subject.ID72 Subject.ID73 Subject.ID74 Subject.ID75 Subject.ID76
## StdDev: 0.1471092 0.1471092 0.1471092 0.1471092 0.1471092
## Subject.ID77 Subject.ID78 Subject.ID79 Subject.ID80 Subject.ID81
## StdDev: 0.1471092 0.1471092 0.1471092 0.1471092 0.1471092
## Subject.ID82 Subject.ID83 Subject.ID84 Subject.ID85 Subject.ID86
## StdDev: 0.1471092 0.1471092 0.1471092 0.1471092 0.1471092
## Subject.ID87 Subject.ID88 Subject.ID89 Subject.ID90 Subject.ID91
## StdDev: 0.1471092 0.1471092 0.1471092 0.1471092 0.1471092
## Subject.ID92 Subject.ID93 Subject.ID94 Subject.ID95 Subject.ID96
## StdDev: 0.1471092 0.1471092 0.1471092 0.1471092 0.1471092
## Subject.ID97 Subject.ID98 Subject.ID99 Subject.ID100 Subject.ID101
## StdDev: 0.1471092 0.1471092 0.1471092 0.1471092 0.1471092
## Subject.ID102 Subject.ID103 Subject.ID104 Subject.ID105 Subject.ID106
## StdDev: 0.1471092 0.1471092 0.1471092 0.1471092 0.1471092
## Subject.ID107 Subject.ID108 Subject.ID109 Subject.ID110 Subject.ID111
## StdDev: 0.1471092 0.1471092 0.1471092 0.1471092 0.1471092
## Subject.ID112 Subject.ID113 Subject.ID114 Subject.ID115 Subject.ID116
## StdDev: 0.1471092 0.1471092 0.1471092 0.1471092 0.1471092
## Subject.ID117 Subject.ID118 Subject.ID119 Subject.ID120 Subject.ID121
## StdDev: 0.1471092 0.1471092 0.1471092 0.1471092 0.1471092
## Subject.ID122 Subject.ID123 Subject.ID124 Subject.ID125 Subject.ID126
## StdDev: 0.1471092 0.1471092 0.1471092 0.1471092 0.1471092
## Subject.ID127 Subject.ID128 Subject.ID129 Subject.ID130 Subject.ID131
## StdDev: 0.1471092 0.1471092 0.1471092 0.1471092 0.1471092
## Subject.ID132 Subject.ID133 Subject.ID134 Subject.ID135 Subject.ID136
## StdDev: 0.1471092 0.1471092 0.1471092 0.1471092 0.1471092
## Subject.ID137 Subject.ID138 Subject.ID139 Subject.ID140 Subject.ID141
## StdDev: 0.1471092 0.1471092 0.1471092 0.1471092 0.1471092
## Subject.ID142 Subject.ID143 Subject.ID144 Subject.ID145 Subject.ID146
## StdDev: 0.1471092 0.1471092 0.1471092 0.1471092 0.1471092
## Subject.ID147 Subject.ID148 Subject.ID149 Subject.ID150 Subject.ID151
## StdDev: 0.1471092 0.1471092 0.1471092 0.1471092 0.1471092
## Subject.ID152 Subject.ID153 Subject.ID154 Subject.ID155 Subject.ID156
## StdDev: 0.1471092 0.1471092 0.1471092 0.1471092 0.1471092
## Subject.ID157 Subject.ID158 Subject.ID159 Subject.ID160 Subject.ID161
## StdDev: 0.1471092 0.1471092 0.1471092 0.1471092 0.1471092
## Subject.ID162 Subject.ID163 Subject.ID164 Subject.ID165 Subject.ID166
## StdDev: 0.1471092 0.1471092 0.1471092 0.1471092 0.1471092
## Subject.ID167 Subject.ID168 Subject.ID169 Subject.ID170 Subject.ID171
## StdDev: 0.1471092 0.1471092 0.1471092 0.1471092 0.1471092
## Subject.ID172 Subject.ID173 Subject.ID174 Subject.ID175 Subject.ID176
## StdDev: 0.1471092 0.1471092 0.1471092 0.1471092 0.1471092
## Subject.ID177 Subject.ID178 Subject.ID179 Subject.ID180 Subject.ID181
## StdDev: 0.1471092 0.1471092 0.1471092 0.1471092 0.1471092
## Subject.ID182 Subject.ID183 Subject.ID184 Subject.ID185 Subject.ID186
## StdDev: 0.1471092 0.1471092 0.1471092 0.1471092 0.1471092
## Subject.ID187 Subject.ID188 Subject.ID189 Subject.ID190 Subject.ID191
## StdDev: 0.1471092 0.1471092 0.1471092 0.1471092 0.1471092
## Subject.ID192 Subject.ID193 Subject.ID194 Subject.ID195 Subject.ID196
## StdDev: 0.1471092 0.1471092 0.1471092 0.1471092 0.1471092
## Subject.ID197 Subject.ID198 Subject.ID199 Subject.ID200 Subject.ID201
## StdDev: 0.1471092 0.1471092 0.1471092 0.1471092 0.1471092
## Subject.ID202 Subject.ID203 Subject.ID204 Subject.ID205 Subject.ID206
## StdDev: 0.1471092 0.1471092 0.1471092 0.1471092 0.1471092
## Subject.ID207 Subject.ID208 Subject.ID209 Subject.ID210 Subject.ID211
## StdDev: 0.1471092 0.1471092 0.1471092 0.1471092 0.1471092
## Subject.ID212 Subject.ID213 Subject.ID214 Subject.ID215 Subject.ID216
## StdDev: 0.1471092 0.1471092 0.1471092 0.1471092 0.1471092
## Subject.ID217 Subject.ID218 Subject.ID219 Subject.ID220 Subject.ID221
## StdDev: 0.1471092 0.1471092 0.1471092 0.1471092 0.1471092
## Subject.ID222 Subject.ID223 Subject.ID224 Subject.ID225 Subject.ID227
## StdDev: 0.1471092 0.1471092 0.1471092 0.1471092 0.1471092
## Subject.ID228 Subject.ID229 Subject.ID230
## StdDev: 0.1471092 0.1471092 0.1471092
##
## Block 4: Degree.Day15, Degree.Day30, Degree.Day45, Degree.Day60, Degree.Day75, Degree.Day90, Degree.Day105, Degree.Day120, Degree.Day135, Degree.Day150, Degree.Day165, Degree.Day180, Degree.Day195, Degree.Day210, Degree.Day285, Degree.Day304, Degree.Day323, Degree.Day342, Degree.Day361, Degree.Day380, Degree.Day399, Degree.Day418, Degree.Day437, Degree.Day456, Degree.Day475, Degree.Day494, Degree.Day513, Degree.Day532, Degree.Day551, Degree.Day570
## Formula: ~0 + Degree.Day | Dummy
## Structure: Multiple of an Identity
## Degree.Day15 Degree.Day30 Degree.Day45 Degree.Day60 Degree.Day75
## StdDev: 0.06349814 0.06349814 0.06349814 0.06349814 0.06349814
## Degree.Day90 Degree.Day105 Degree.Day120 Degree.Day135 Degree.Day150
## StdDev: 0.06349814 0.06349814 0.06349814 0.06349814 0.06349814
## Degree.Day165 Degree.Day180 Degree.Day195 Degree.Day210 Degree.Day285
## StdDev: 0.06349814 0.06349814 0.06349814 0.06349814 0.06349814
## Degree.Day304 Degree.Day323 Degree.Day342 Degree.Day361 Degree.Day380
## StdDev: 0.06349814 0.06349814 0.06349814 0.06349814 0.06349814
## Degree.Day399 Degree.Day418 Degree.Day437 Degree.Day456 Degree.Day475
## StdDev: 0.06349814 0.06349814 0.06349814 0.06349814 0.06349814
## Degree.Day494 Degree.Day513 Degree.Day532 Degree.Day551 Degree.Day570
## StdDev: 0.06349814 0.06349814 0.06349814 0.06349814 0.06349814
## Residual
## StdDev: 0.2026764
##
## Correlation Structure: AR(1)
## Formula: ~1 | Dummy
## Parameter estimate(s):
## Phi
## 0.3510209
## Fixed effects: MO2.NMS ~ Mass + bSpline * Offspring.Acc * Father.Acc * Mother.Acc
## Value Std.Error DF
## (Intercept) -0.8447864 0.11289296 2554
## Mass 0.2049529 0.01670045 2554
## bSpline 1.1932080 0.09159720 2554
## Offspring.Acc2 0.0173138 0.11502032 2554
## Father.Acc2 0.4360736 0.14145323 2554
## Mother.Acc2 0.2226406 0.10883690 2554
## bSpline:Offspring.Acc2 -0.0028834 0.12384042 2554
## bSpline:Father.Acc2 -0.3593654 0.13311546 2554
## Offspring.Acc2:Father.Acc2 -0.0476710 0.15979645 2554
## bSpline:Mother.Acc2 -0.1737977 0.12144948 2554
## Offspring.Acc2:Mother.Acc2 0.0112290 0.15280483 2554
## Father.Acc2:Mother.Acc2 -0.2187508 0.15867081 2554
## bSpline:Offspring.Acc2:Father.Acc2 -0.0482553 0.17574242 2554
## bSpline:Offspring.Acc2:Mother.Acc2 -0.1292361 0.16544425 2554
## bSpline:Father.Acc2:Mother.Acc2 0.1691838 0.18017240 2554
## Offspring.Acc2:Father.Acc2:Mother.Acc2 -0.1522449 0.21895744 2554
## bSpline:Offspring.Acc2:Father.Acc2:Mother.Acc2 0.3038288 0.23939520 2554
## t-value p-value
## (Intercept) -7.483075 0.0000
## Mass 12.272296 0.0000
## bSpline 13.026686 0.0000
## Offspring.Acc2 0.150528 0.8804
## Father.Acc2 3.082811 0.0021
## Mother.Acc2 2.045635 0.0409
## bSpline:Offspring.Acc2 -0.023283 0.9814
## bSpline:Father.Acc2 -2.699652 0.0070
## Offspring.Acc2:Father.Acc2 -0.298323 0.7655
## bSpline:Mother.Acc2 -1.431029 0.1525
## Offspring.Acc2:Mother.Acc2 0.073486 0.9414
## Father.Acc2:Mother.Acc2 -1.378646 0.1681
## bSpline:Offspring.Acc2:Father.Acc2 -0.274580 0.7837
## bSpline:Offspring.Acc2:Mother.Acc2 -0.781146 0.4348
## bSpline:Father.Acc2:Mother.Acc2 0.939010 0.3478
## Offspring.Acc2:Father.Acc2:Mother.Acc2 -0.695317 0.4869
## bSpline:Offspring.Acc2:Father.Acc2:Mother.Acc2 1.269152 0.2045
## Correlation:
## (Intr) Mass bSplin Off.A2
## Mass -0.475
## bSpline -0.618 -0.002
## Offspring.Acc2 -0.487 -0.081 0.607
## Father.Acc2 -0.647 0.095 0.494 0.397
## Mother.Acc2 -0.559 0.029 0.641 0.532
## bSpline:Offspring.Acc2 0.460 -0.006 -0.740 -0.841
## bSpline:Father.Acc2 0.425 0.005 -0.690 -0.420
## Offspring.Acc2:Father.Acc2 0.365 0.000 -0.438 -0.694
## bSpline:Mother.Acc2 0.466 0.001 -0.754 -0.458
## Offspring.Acc2:Mother.Acc2 0.382 0.011 -0.456 -0.735
## Father.Acc2:Mother.Acc2 0.367 -0.006 -0.441 -0.357
## bSpline:Offspring.Acc2:Father.Acc2 -0.325 0.002 0.522 0.596
## bSpline:Offspring.Acc2:Mother.Acc2 -0.347 0.009 0.554 0.629
## bSpline:Father.Acc2:Mother.Acc2 -0.314 -0.004 0.510 0.310
## Offspring.Acc2:Father.Acc2:Mother.Acc2 -0.252 -0.027 0.321 0.507
## bSpline:Offspring.Acc2:Father.Acc2:Mother.Acc2 0.240 -0.004 -0.385 -0.437
## Fth.A2 Mth.A2 bSp:O.A2 bSp:F.A2
## Mass
## bSpline
## Offspring.Acc2
## Father.Acc2
## Mother.Acc2 0.425
## bSpline:Offspring.Acc2 -0.366 -0.474
## bSpline:Father.Acc2 -0.714 -0.442 0.510
## Offspring.Acc2:Father.Acc2 -0.611 -0.373 0.607 0.632
## bSpline:Mother.Acc2 -0.372 -0.846 0.558 0.519
## Offspring.Acc2:Mother.Acc2 -0.298 -0.710 0.633 0.315
## Father.Acc2:Mother.Acc2 -0.608 -0.672 0.326 0.636
## bSpline:Offspring.Acc2:Father.Acc2 0.541 0.336 -0.707 -0.757
## bSpline:Offspring.Acc2:Mother.Acc2 0.274 0.621 -0.749 -0.381
## bSpline:Father.Acc2:Mother.Acc2 0.527 0.572 -0.377 -0.738
## Offspring.Acc2:Father.Acc2:Mother.Acc2 0.438 0.486 -0.444 -0.462
## bSpline:Offspring.Acc2:Father.Acc2:Mother.Acc2 -0.397 -0.431 0.519 0.556
## Of.A2:F.A2 bS:M.A O.A2:M F.A2:M
## Mass
## bSpline
## Offspring.Acc2
## Father.Acc2
## Mother.Acc2
## bSpline:Offspring.Acc2
## bSpline:Father.Acc2
## Offspring.Acc2:Father.Acc2
## bSpline:Mother.Acc2 0.330
## Offspring.Acc2:Mother.Acc2 0.517 0.602
## Father.Acc2:Mother.Acc2 0.537 0.581 0.478
## bSpline:Offspring.Acc2:Father.Acc2 -0.857 -0.393 -0.448 -0.482
## bSpline:Offspring.Acc2:Mother.Acc2 -0.453 -0.734 -0.844 -0.426
## bSpline:Father.Acc2:Mother.Acc2 -0.467 -0.675 -0.407 -0.859
## Offspring.Acc2:Father.Acc2:Mother.Acc2 -0.727 -0.422 -0.692 -0.725
## bSpline:Offspring.Acc2:Father.Acc2:Mother.Acc2 0.629 0.509 0.586 0.647
## bSp:O.A2:F.A2 bS:O.A2:M bS:F.A2:
## Mass
## bSpline
## Offspring.Acc2
## Father.Acc2
## Mother.Acc2
## bSpline:Offspring.Acc2
## bSpline:Father.Acc2
## Offspring.Acc2:Father.Acc2
## bSpline:Mother.Acc2
## Offspring.Acc2:Mother.Acc2
## Father.Acc2:Mother.Acc2
## bSpline:Offspring.Acc2:Father.Acc2
## bSpline:Offspring.Acc2:Mother.Acc2 0.528
## bSpline:Father.Acc2:Mother.Acc2 0.559 0.495
## Offspring.Acc2:Father.Acc2:Mother.Acc2 0.626 0.590 0.623
## bSpline:Offspring.Acc2:Father.Acc2:Mother.Acc2 -0.734 -0.693 -0.753
## O.A2:F.A2:
## Mass
## bSpline
## Offspring.Acc2
## Father.Acc2
## Mother.Acc2
## bSpline:Offspring.Acc2
## bSpline:Father.Acc2
## Offspring.Acc2:Father.Acc2
## bSpline:Mother.Acc2
## Offspring.Acc2:Mother.Acc2
## Father.Acc2:Mother.Acc2
## bSpline:Offspring.Acc2:Father.Acc2
## bSpline:Offspring.Acc2:Mother.Acc2
## bSpline:Father.Acc2:Mother.Acc2
## Offspring.Acc2:Father.Acc2:Mother.Acc2
## bSpline:Offspring.Acc2:Father.Acc2:Mother.Acc2 -0.849
##
## Standardized Within-Group Residuals:
## Min Q1 Med Q3 Max
## -3.02357916 -0.61273768 -0.06675436 0.56340609 3.77352550
##
## Number of Observations: 2571
## Number of Groups: 1
# Model output with corrected degrees of freedom and p-values
print(coef_output)
## Value Std.Error
## (Intercept) -0.844786431 0.11289296
## Mass 0.204952884 0.01670045
## bSpline 1.193208038 0.09159720
## Offspring.Acc2 0.017313773 0.11502032
## Father.Acc2 0.436073634 0.14145323
## Mother.Acc2 0.222640621 0.10883690
## bSpline:Offspring.Acc2 -0.002883363 0.12384042
## bSpline:Father.Acc2 -0.359365393 0.13311546
## Offspring.Acc2:Father.Acc2 -0.047671022 0.15979645
## bSpline:Mother.Acc2 -0.173797699 0.12144948
## Offspring.Acc2:Mother.Acc2 0.011229015 0.15280483
## Father.Acc2:Mother.Acc2 -0.218750846 0.15867081
## bSpline:Offspring.Acc2:Father.Acc2 -0.048255289 0.17574242
## bSpline:Offspring.Acc2:Mother.Acc2 -0.129236061 0.16544425
## bSpline:Father.Acc2:Mother.Acc2 0.169183759 0.18017240
## Offspring.Acc2:Father.Acc2:Mother.Acc2 -0.152244934 0.21895744
## bSpline:Offspring.Acc2:Father.Acc2:Mother.Acc2 0.303828827 0.23939520
## DF t.value
## (Intercept) 19.04531 -7.48307473
## Mass 160.99689 12.27229627
## bSpline 2338.65930 13.02668634
## Offspring.Acc2 415.63258 0.15052795
## Father.Acc2 11.41576 3.08281137
## Mother.Acc2 213.08059 2.04563541
## bSpline:Offspring.Acc2 2340.43281 -0.02328290
## bSpline:Father.Acc2 2338.81306 -2.69965173
## Offspring.Acc2:Father.Acc2 767.14285 -0.29832341
## bSpline:Mother.Acc2 2337.24613 -1.43102876
## Offspring.Acc2:Mother.Acc2 470.20003 0.07348599
## Father.Acc2:Mother.Acc2 834.12226 -1.37864582
## bSpline:Offspring.Acc2:Father.Acc2 2343.31538 -0.27457964
## bSpline:Offspring.Acc2:Mother.Acc2 2345.15620 -0.78114569
## bSpline:Father.Acc2:Mother.Acc2 2338.46722 0.93901039
## Offspring.Acc2:Father.Acc2:Mother.Acc2 557.56944 -0.69531748
## bSpline:Offspring.Acc2:Father.Acc2:Mother.Acc2 2344.32382 1.26915172
## p.value
## (Intercept) 4.387906e-07
## Mass 7.360375e-25
## bSpline 1.684939e-37
## Offspring.Acc2 8.804211e-01
## Father.Acc2 1.000747e-02
## Mother.Acc2 4.202078e-02
## bSpline:Offspring.Acc2 9.814266e-01
## bSpline:Father.Acc2 6.991185e-03
## Offspring.Acc2:Father.Acc2 7.655371e-01
## bSpline:Mother.Acc2 1.525557e-01
## Offspring.Acc2:Mother.Acc2 9.414506e-01
## Father.Acc2:Mother.Acc2 1.683736e-01
## bSpline:Offspring.Acc2:Father.Acc2 7.836634e-01
## bSpline:Offspring.Acc2:Mother.Acc2 4.347957e-01
## bSpline:Father.Acc2:Mother.Acc2 3.478224e-01
## Offspring.Acc2:Father.Acc2:Mother.Acc2 4.871460e-01
## bSpline:Offspring.Acc2:Father.Acc2:Mother.Acc2 2.045129e-01
summary(ModelTMB.AR)
## Family: gaussian ( identity )
## Formula:
## MO2.NMS ~ Mass + bSpline * Offspring.Acc * Father.Acc * Mother.Acc +
## (1 | Father.ID) + (1 | Mother.ID) + (1 | Subject.ID) + (1 |
## Degree.Day) + ar1(Time + 0 | AR_Group)
## Data: BTChallenge_Clean
##
## AIC BIC logLik deviance df.resid
## -993.6 -853.2 520.8 -1041.6 2542
##
## Random effects:
##
## Conditional model:
## Groups Name Variance Std.Dev. Corr
## Father.ID (Intercept) 1.977e-03 4.446e-02
## Mother.ID (Intercept) 1.530e-13 3.911e-07
## Subject.ID (Intercept) 3.160e-10 1.778e-05
## Degree.Day (Intercept) 2.000e-03 4.472e-02
## AR_Group Time1 4.335e-02 2.082e-01 0.89 (ar1)
## Residual 1.998e-02 1.414e-01
## Number of obs: 2566, groups:
## Father.ID, 4; Mother.ID, 4; Subject.ID, 229; Degree.Day, 30; AR_Group, 229
##
## Dispersion estimate for gaussian family (sigma^2): 0.02
##
## Conditional model:
## Estimate Std. Error z value
## (Intercept) -0.78843 0.11501 -6.855
## Mass 0.20446 0.01690 12.096
## bSpline 1.12275 0.10919 10.283
## Offspring.Acc2 -0.03014 0.13249 -0.227
## Father.Acc2 0.37947 0.14590 2.601
## Mother.Acc2 0.11575 0.12572 0.921
## bSpline:Offspring.Acc2 0.03452 0.14663 0.235
## bSpline:Father.Acc2 -0.28727 0.15767 -1.822
## Offspring.Acc2:Father.Acc2 -0.06369 0.18608 -0.342
## bSpline:Mother.Acc2 -0.01036 0.14333 -0.072
## Offspring.Acc2:Mother.Acc2 0.06761 0.17570 0.385
## Father.Acc2:Mother.Acc2 -0.13267 0.18465 -0.718
## bSpline:Offspring.Acc2:Father.Acc2 -0.01277 0.20730 -0.062
## bSpline:Offspring.Acc2:Mother.Acc2 -0.19074 0.19490 -0.979
## bSpline:Father.Acc2:Mother.Acc2 0.04054 0.21134 0.192
## Offspring.Acc2:Father.Acc2:Mother.Acc2 -0.14540 0.25177 -0.578
## bSpline:Offspring.Acc2:Father.Acc2:Mother.Acc2 0.30242 0.28061 1.078
## Pr(>|z|)
## (Intercept) 7.13e-12 ***
## Mass < 2e-16 ***
## bSpline < 2e-16 ***
## Offspring.Acc2 0.8201
## Father.Acc2 0.0093 **
## Mother.Acc2 0.3572
## bSpline:Offspring.Acc2 0.8139
## bSpline:Father.Acc2 0.0685 .
## Offspring.Acc2:Father.Acc2 0.7322
## bSpline:Mother.Acc2 0.9424
## Offspring.Acc2:Mother.Acc2 0.7004
## Father.Acc2:Mother.Acc2 0.4725
## bSpline:Offspring.Acc2:Father.Acc2 0.9509
## bSpline:Offspring.Acc2:Mother.Acc2 0.3277
## bSpline:Father.Acc2:Mother.Acc2 0.8479
## Offspring.Acc2:Father.Acc2:Mother.Acc2 0.5636
## bSpline:Offspring.Acc2:Father.Acc2:Mother.Acc2 0.2812
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
View(data.frame(
"Term" = rownames(coef_output),
"lme_Beta" = coef_output$Value, "TMB_Beta" = as.numeric(summary(ModelTMB.AR)$coefficient$cond[, 1]),
"lme_p" = coef_output$p.value, "TMB_p" = as.numeric(summary(ModelTMB.AR)$coefficient$cond[, 4])
) %>%
mutate("Diff" = ifelse(lme_p < 0.05 & TMB_p > 0.05, "*",
ifelse(TMB_p < 0.05 & lme_p > 0.05, "*", " ")
)))
# Some subtle differences between models
# Now using Kenward-Roger approximation
require(pbkrtest)
get_kr_df <- function(model_object) {
L <- diag(rep(1, length(fixef(model_object))))
L <- as.data.frame(L)
out <- purrr::map_dbl(L, pbkrtest::get_Lb_ddf, object = model_object)
names(out) <- names(fixef(model_object))
out
}
KR_DF <- get_kr_df(DF.Model)
p_vals <- c()
for (i in 1:length(t_vals)) {
if (t_vals[i] > 0) {
p_vals[i] <- 2 * (pt(t_vals[i], df = as.numeric(KR_DF)[i], lower = FALSE))
} else if (t_vals[i] < 0) {
p_vals[i] <- 2 * (pt(t_vals[i], df = as.numeric(KR_DF)[i], lower = TRUE))
}
}
DF_coef_output <- data.frame(summary(Model.AR)$tTab)
DF_coef_output$p.value <- p_vals
DF_coef_output$DF <- KR_DF
View(data.frame(
"Term" = rownames(coef_output),
"lme_BetaSat" = coef_output$Value, "lme_BetaKR" = DF_coef_output$Value,
"TMB_Beta" = as.numeric(summary(ModelTMB.AR)$coefficient$cond[, 1]),
"lme_pSat" = coef_output$p.value, "lme_pKR" = DF_coef_output$p.value,
"TMB_p" = as.numeric(summary(ModelTMB.AR)$coefficient$cond[, 4])
) %>%
mutate("Diff" = ifelse(lme_pKR < 0.05 & TMB_p > 0.05, "*",
ifelse(TMB_p < 0.05 & lme_pKR > 0.05, "*", " ")
)))
# Satherwaitte and KR outputs are very similar - probably owing to such large DF estimates.
# Estimating dfs according to sample sizes of offspring alone.
p_vals <- c()
for (i in 1:length(t_vals)) {
if (t_vals[i] > 0) {
p_vals[i] <- 2 * (pt(t_vals[i], df = 228, lower = FALSE))
} else if (t_vals[i] < 0) {
p_vals[i] <- 2 * (pt(t_vals[i], df = 228, lower = TRUE))
}
}
Simple_coef_output <- data.frame(summary(Model.AR)$tTab)
Simple_coef_output$p.value <- p_vals
Simple_coef_output$DF <- 228
View(data.frame(
"Term" = rownames(coef_output),
"lme_Beta" = coef_output$Value,
"TMB_Beta" = as.numeric(summary(ModelTMB.AR)$coefficient$cond[, 1]),
"lme_pSat" = coef_output$p.value, "lme_pKR" = DF_coef_output$p.value,
"lme_pSimple" = Simple_coef_output$p.value,
"TMB_p" = as.numeric(summary(ModelTMB.AR)$coefficient$cond[, 4])
) %>%
mutate("Diff" = ifelse(lme_pSimple < 0.05 & lme_pKR > 0.05, "*",
ifelse(lme_pKR < 0.05 & lme_pSimple > 0.05, "*", " ")
)))
# Printing off data frame
{
T1 <- coef_output[, c("Value", "DF", "t.value", "p.value")] %>%
mutate("Method" = "Satterthwaite", "Term" = rownames(coef_output)) %>%
rename(
"Test_Statistic" = "t.value", "Estimate" = "Value",
"p" = "p.value"
) %>%
select(Term, Estimate, Method, DF, Test_Statistic, p)
T2 <- DF_coef_output[, c("Value", "DF", "t.value", "p.value")] %>%
mutate("Method" = "Kenward-Roger", "Term" = rownames(coef_output)) %>%
rename(
"Test_Statistic" = "t.value", "Estimate" = "Value",
"p" = "p.value"
) %>%
select(Term, Estimate, Method, DF, Test_Statistic, p)
T3 <- as.data.frame(summary(ModelTMB.AR)$coefficient$cond)[, c("Estimate", "z value", "Pr(>|z|)")] %>%
mutate(
"Method" = "z-value", "DF" = NA,
"Term" = rownames(as.data.frame(summary(ModelTMB.AR)$coefficient$cond))
) %>%
rename(
"Test_Statistic" = "z value",
"p" = "Pr(>|z|)"
) %>%
select(Term, Estimate, Method, DF, Test_Statistic, p)
T4 <- Simple_coef_output[, c("Value", "DF", "t.value", "p.value")] %>%
mutate("Method" = "n-1", "Term" = rownames(coef_output)) %>%
rename(
"Test_Statistic" = "t.value", "Estimate" = "Value",
"p" = "p.value"
) %>%
select(Term, Estimate, Method, DF, Test_Statistic, p)
}
All <- rbind(T1, T2, T3, T4) %>%
mutate(Method = factor(Method,
levels = c("Satterthwaite", "Kenward-Roger", "n-1", "z-value")
)) %>%
arrange(Term, Method) %>%
mutate("p < alpha" = ifelse(p < 0.05, "*", ""))
write.csv(All, "Model_Comparison_Final.csv", row.names = FALSE)
Calculating marginal and conditional r-squared values
r.squaredGLMM(Model.AR)
## R2m R2c
## [1,] 0.3932383 0.6577368
# Quite low in comparison to models in Penney et al (In Press), but somewhat reasonable.
Producing marginal mean plots for lme model
require("emmeans")
{
EM_Means_Off <- emmeans::emmip(Model.AR, ~ bSpline | Offspring.Acc, cov.reduce = FALSE, nesting = NULL, nesting.order = FALSE)
EM_Means_Mat <- emmeans::emmip(Model.AR, ~ bSpline | Mother.Acc, cov.reduce = FALSE)
EM_Means_Pat <- emmeans::emmip(Model.AR, ~ bSpline | Father.Acc ~ bSpline, cov.reduce = FALSE)
EM_Means_Off_Pat <- emmeans::emmip(Model.AR, Offspring.Acc ~ bSpline | Father.Acc, cov.reduce = FALSE)
EM_Means_Off_Mat <- emmeans::emmip(Model.AR, Offspring.Acc ~ bSpline | Mother.Acc, cov.reduce = FALSE)
EM_Means_Mass <- emmeans::emmip(Model.AR, ~Mass, at = list(Mass = seq(min(Dat$Mass, na.rm = T), max(Dat$Mass, na.rm = T), 0.1)), CIs = TRUE, type = "response", data = Dat)
}
Raw <- as.data.frame(emmeans::emmeans(Model.AR, specs = "bSpline", by = "Offspring.Acc", cov.reduce = FALSE))
bSpline_vals <- Raw$bSpline
Challenge_Temp <- c()
for (i in 1:length(bSpline_vals)) {
Challenge_Temp[i] <- BTChallenge[c(which(BTChallenge$bSpline == bSpline_vals[i])), 9]
}
# Removing replicates
Challenge_Temp <- unique(Challenge_Temp)
EM_Off <- as.data.frame(EM_Means_Off$data)
EM_Off$Challenge <- sort(rep(Challenge_Temp, 2))
# Removing prodictions for warm acclimated offspring at cold temperatures
to_rm <- c(which(EM_Off$Offspring.Acc == "2" & EM_Off$Challenge < 19))
EM_Off <- EM_Off[-to_rm, ]
EM_Off$Upper.CI <- EM_Off$yvar + EM_Off$SE * 1.96
EM_Off$Lower.CI <- EM_Off$yvar - EM_Off$SE * 1.96
EM_Plot_Off_Ribbon <- ggplot(EM_Off, aes(x = Challenge, y = yvar, fill = Offspring.Acc, linetype = Offspring.Acc)) +
theme_bw() +
geom_smooth(method = "gam", formula = y ~ s(x, k = 7, bs = "cs"), se = FALSE, colour = "black") +
geom_ribbon(mapping = aes(x = Challenge, ymin = c(Lower.CI), ymax = c(Upper.CI)), alpha = 0.3, colour = NA) +
scale_fill_manual(values = c("dodgerblue3", "magenta"), labels = c("Cold", "Warm")) +
scale_linetype_manual(values = c("solid", "dashed"), labels = c("Cold", "Warm")) +
labs(fill = "Offspring Acclimation \n Temperature", linetype = "Offspring Acclimation \n Temperature") +
theme(axis.text = element_text(size = 8, colour = "black", family = "Noto Sans")) +
labs(x = "Challenge Temperature (°C)", y = expression("Resting MO"[2] * " (mg O"[2] * "/h)")) +
theme(panel.grid.major = element_blank()) +
xlim(c(15, 30))
EM_Mat <- as.data.frame(EM_Means_Mat$data)
EM_Mat$Challenge <- sort(rep(Challenge_Temp, 2))
EM_Mat$Upper.CI <- EM_Mat$yvar + EM_Mat$SE * 1.96
EM_Mat$Lower.CI <- EM_Mat$yvar - EM_Mat$SE * 1.96
EM_Plot_Mat_Ribbon <- ggplot(EM_Mat, aes(x = Challenge, y = yvar, linetype = Mother.Acc, fill = Mother.Acc)) +
theme_bw() +
geom_smooth(method = "gam", formula = y ~ s(x, k = 7, bs = "cs"), se = FALSE, colour = "black") +
geom_ribbon(mapping = aes(x = Challenge, ymin = c(Lower.CI), ymax = c(Upper.CI)), alpha = 0.3, colour = NA) +
scale_fill_manual(values = c("dodgerblue3", "magenta"), labels = c("Cold", "Warm")) +
scale_linetype_manual(values = c("solid", "dashed"), labels = c("Cold", "Warm")) +
labs(fill = "Mother Acclimation \n Temperature", linetype = "Mother Acclimation \n Temperature") +
theme(axis.text = element_text(size = 8, colour = "black", family = "Noto Sans")) +
labs(x = "Challenge Temperature (°C)", y = expression("Resting MO"[2] * " (mg O"[2] * "/h)")) +
theme(panel.grid.major = element_blank()) +
xlim(c(15, 30))
EM_Pat <- as.data.frame(EM_Means_Pat$data)
EM_Pat$Challenge <- sort(rep(Challenge_Temp, 2))
EM_Pat$Upper.CI <- EM_Pat$yvar + EM_Pat$SE * 1.96
EM_Pat$Lower.CI <- EM_Pat$yvar - EM_Pat$SE * 1.96
EM_Plot_Pat_Ribbon <- ggplot(EM_Pat, aes(x = Challenge, y = yvar, fill = Father.Acc, linetype = Father.Acc)) +
theme_bw() +
geom_smooth(method = "gam", formula = y ~ s(x, k = 7, bs = "cs"), se = FALSE, colour = "black") +
geom_ribbon(mapping = aes(x = Challenge, ymin = c(Lower.CI), ymax = c(Upper.CI)), alpha = 0.3, colour = NA) +
scale_fill_manual(values = c("dodgerblue3", "magenta"), labels = c("Cold", "Warm")) +
scale_linetype_manual(values = c("solid", "dashed"), labels = c("Cold", "Warm")) +
labs(fill = "Father Acclimation \n Temperature", linetype = "Father Acclimation \n Temperature") +
theme(axis.text = element_text(size = 8, colour = "black", family = "Noto Sans")) +
labs(x = "Challenge Temperature (°C)", y = expression("Resting MO"[2] * " (mg O"[2] * "/h)")) +
theme(panel.grid.major = element_blank()) +
xlim(c(15, 30))
Off_Pat <- as.data.frame(EM_Means_Off_Pat$data)
Off_Pat$Challenge <- sort(rep(Challenge_Temp, 4))
Off_Pat$Upper.CI <- Off_Pat$yvar + Off_Pat$SE * 1.96
Off_Pat$Lower.CI <- Off_Pat$yvar - Off_Pat$SE * 1.96
to_rm <- c(which(Off_Pat$Offspring.Acc == "2" & Off_Pat$Challenge < 19))
Off_Pat <- Off_Pat[-to_rm, ]
facet_names <- c(`1` = "Cold Fathers", `2` = "Warm Fathers")
facet_names_new <- c(`1` = "Cold Offspring", `2` = "Warm Offspring")
EM_Plot_Off_Pat_Ribbon <- ggplot(Off_Pat, aes(x = Challenge, y = yvar, fill = Father.Acc, linetype = Father.Acc)) +
theme_bw() +
geom_smooth(method = "gam", formula = y ~ s(x, k = 7, bs = "cs"), se = FALSE, colour = "black") +
geom_ribbon(mapping = aes(x = Challenge, ymin = c(Lower.CI), ymax = c(Upper.CI)), alpha = 0.3, colour = NA) +
scale_fill_manual(values = c("dodgerblue3", "magenta"), labels = c("Cold", "Warm")) +
scale_linetype_manual(values = c("solid", "dashed"), labels = c("Cold", "Warm")) +
labs(linetype = "Father Acclimation \n Temperature", fill = "Father Acclimation \n Temperature") +
theme(axis.text = element_text(size = 8, colour = "black", family = "Noto Sans")) +
facet_wrap(~Offspring.Acc, labeller = as_labeller(facet_names_new)) +
labs(x = "Challenge Temperature (°C)", y = expression("Resting MO"[2] * " (mg O"[2] * "/h)")) +
theme(strip.text.x = element_text(family = "Noto Sans", colour = "black"), panel.grid.major = element_blank()) +
xlim(c(15, 30))
Off_Mat <- as.data.frame(EM_Means_Off_Mat$data)
Off_Mat$Challenge <- sort(rep(Challenge_Temp, 4))
Off_Mat$Upper.CI <- Off_Mat$yvar + Off_Mat$SE * 1.96
Off_Mat$Lower.CI <- Off_Mat$yvar - Off_Mat$SE * 1.96
to_rm <- c(which(Off_Mat$Offspring.Acc == "2" & Off_Mat$Challenge < 19))
Off_Mat <- Off_Mat[-to_rm, ]
facet_names <- c(`1` = "Cold Mothers", `2` = "Warm Mothers")
EM_Plot_Off_Mat_Ribbon <- ggplot(Off_Mat, aes(x = Challenge, y = yvar, fill = Mother.Acc, linetype = Mother.Acc)) +
theme_bw() +
geom_smooth(method = "gam", formula = y ~ s(x, k = 7, bs = "cs"), se = FALSE, colour = "black") +
geom_ribbon(mapping = aes(x = Challenge, ymin = c(Lower.CI), ymax = c(Upper.CI)), alpha = 0.3, colour = NA) +
scale_fill_manual(values = c("dodgerblue3", "magenta"), labels = c("Cold", "Warm")) +
scale_linetype_manual(values = c("solid", "dashed"), labels = c("Cold", "Warm")) +
labs(linetype = "Mother Acclimation \n Temperature", fill = "Mother Acclimation \n Temperature") +
theme(axis.text = element_text(size = 8, colour = "black", family = "Noto Sans")) +
facet_wrap(~Offspring.Acc, labeller = as_labeller(facet_names_new)) +
labs(x = "Challenge Temperature (°C)", y = expression("Resting MO"[2] * " (mg O"[2] * "/h)")) +
theme(strip.text.x = element_text(family = "Noto Sans", colour = "black"), panel.grid.major = element_blank()) +
xlim(c(15, 30))
EM_Mass <- as.data.frame(EM_Means_Mass$data)
EM_Mass$Upper.CI <- EM_Mass$yvar + EM_Mass$SE * 1.96
EM_Mass$Lower.CI <- EM_Mass$yvar - EM_Mass$SE * 1.96
EM_Plot_Mass <- ggplot(EM_Mass, aes(x = Mass, y = yvar)) +
theme_bw() +
geom_smooth(method = "gam", formula = y ~ s(x, k = 7, bs = "cs"), se = FALSE, colour = "black") +
geom_ribbon(mapping = aes(x = Mass, ymin = c(Lower.CI), ymax = c(Upper.CI)), alpha = 0.3, colour = NA) +
xlab("Mass (g)") +
theme(axis.text = element_text(size = 8, colour = "black", family = "Noto Sans")) +
labs(x = "Mass (g)", y = expression("MO"[2] * " (mg O"[2] * "/h)")) +
theme(panel.grid.major = element_blank())
grid.arrange(EM_Plot_Off_Ribbon, EM_Plot_Mat_Ribbon, EM_Plot_Pat_Ribbon, EM_Plot_Off_Mat_Ribbon, EM_Plot_Off_Pat_Ribbon, EM_Plot_Mass, nrow = 3)
Spline_plot <- BTChallenge %>%
ggplot(aes(x = Challenge.Temp, y = MO2.NMS)) +
geom_point(pch = 21, size = 3, alpha = 0.3, colour = "black", fill = "royalblue3") +
geom_smooth(
mapping = aes(x = Challenge.Temp, y = bSpline),
colour = "black", linetype = "dashed"
) +
labs(x = "Challenge Temperature (°C)", y = expression("MO"[2] * " (mg O"[2] * "/h)")) +
theme_bw() +
theme(axis.text = element_text(size = 8, colour = "black", family = "Noto Sans"))
print(Spline_plot)
# Producing plots to save
dev.new()
EM_Plot_Off_Ribbon
showtext_opts(dpi = 800)
showtext_auto()
ggsave("/home/joshk/Documents/Trout/Brook Trout Challenge Data/Offspring_alone.jpeg", EM_Plot_Off_Ribbon, dpi = 800, width = 9.45, height = 6.72)
showtext_auto(enable = FALSE)
dev.new()
EM_Plot_Mat_Ribbon
showtext_opts(dpi = 800)
showtext_auto()
ggsave("/home/joshk/Documents/Trout/Brook Trout Challenge Data/Mothers_alone.jpeg", EM_Plot_Mat_Ribbon, dpi = 800, width = 9.45, height = 6.72)
showtext_auto(enable = FALSE)
dev.new()
EM_Plot_Pat_Ribbon
showtext_opts(dpi = 800)
showtext_auto()
ggsave("/home/joshk/Documents/Trout/Brook Trout Challenge Data/Fathers_alone.jpeg", EM_Plot_Pat_Ribbon, dpi = 800, width = 9.45, height = 6.72)
showtext_auto(enable = FALSE)
dev.new()
EM_Plot_Off_Mat_Ribbon
showtext_opts(dpi = 800)
showtext_auto()
ggsave("/home/joshk/Documents/Trout/Brook Trout Challenge Data/Offspring_by_mothers.jpeg", EM_Plot_Off_Mat_Ribbon, dpi = 800, width = 10, height = 6.72)
showtext_auto(enable = FALSE)
dev.new()
EM_Plot_Off_Pat_Ribbon
showtext_opts(dpi = 800)
showtext_auto()
ggsave("/home/joshk/Documents/Trout/Brook Trout Challenge Data/Offspring_by_fathers.jpeg", EM_Plot_Off_Pat_Ribbon, dpi = 800, width = 10, height = 6.72)
showtext_auto(enable = FALSE)
dev.new()
EM_Plot_Mass
showtext_opts(dpi = 800)
showtext_auto()
ggsave("/home/joshk/Documents/Trout/Brook Trout Challenge Data/Mass_alone.jpeg", EM_Plot_Mass, dpi = 800, width = 9.45, height = 6.72)
showtext_auto(enable = FALSE)
dev.new()
Spline_plot
showtext_opts(dpi = 800)
showtext_auto()
ggsave("/home/joshk/Documents/Trout/Brook Trout Challenge Data/Spline_plot.jpeg", Spline_plot, dpi = 800, width = 9.45, height = 6.72)
showtext_auto(enable = FALSE)
## With dots
EM_Plot_Off_Ribbon <- ggplot(EM_Off, aes(x = Challenge, y = yvar, fill = Offspring.Acc, linetype = Offspring.Acc)) +
theme_bw() +
geom_smooth(method = "gam", formula = y ~ s(x, k = 7, bs = "cs"), se = FALSE, colour = "black") +
geom_ribbon(mapping = aes(x = Challenge, ymin = c(Lower.CI), ymax = c(Upper.CI)), alpha = 0.3, colour = NA) +
scale_linetype_manual(values = c("solid", "dashed"), labels = c("Cold", "Warm")) +
scale_fill_manual(values = c("dodgerblue3", "magenta"), labels = c("Cold", "Warm")) +
labs(fill = "Offspring Acclimation \n Temperature", linetype = "Offspring Acclimation \n Temperature") +
theme(axis.text = element_text(size = 8, colour = "black", family = "Noto Sans")) +
labs(x = "Challenge Temperature (°C)", y = expression("Resting MO"[2] * " (mg O"[2] * "/h)")) +
theme(panel.grid.major = element_blank()) +
xlim(c(15, 30)) +
stat_summary(geom = "point", fun = "mean", size = 2, colour = "black", pch = 21, alpha = 0.2, data = Dat, aes(x = Challenge.Temp, y = MO2.NMS))
EM_Plot_Mat_Ribbon <- ggplot(EM_Mat, aes(x = Challenge, y = yvar, colour = Mother.Acc, fill = Mother.Acc)) +
theme_bw() +
geom_smooth(method = "gam", formula = y ~ s(x, k = 7, bs = "cs"), se = FALSE) +
geom_ribbon(mapping = aes(x = Challenge, ymin = c(Lower.CI), ymax = c(Upper.CI)), alpha = 0.3, colour = NA) +
scale_fill_manual(values = c("dodgerblue3", "magenta"), labels = c("Cold", "Warm")) +
scale_colour_manual(values = c("dodgerblue3", "magenta"), labels = c("Cold", "Warm")) +
labs(fill = "Mother Acclimation \n Temperature", colour = "Mother Acclimation \n Temperature") +
theme(axis.text = element_text(size = 8, colour = "black", family = "Noto Sans")) +
labs(x = "Challenge Temperature (°C)", y = expression("Resting MO"[2] * " (mg O"[2] * "/h)")) +
theme(panel.grid.major = element_blank()) +
xlim(c(15, 30)) +
stat_summary(geom = "point", fun = "mean", size = 2, colour = "black", pch = 21, alpha = 0.2, data = Dat, aes(x = Challenge.Temp, y = MO2.NMS))
EM_Plot_Pat_Ribbon <- ggplot(EM_Pat, aes(x = Challenge, y = yvar, fill = Father.Acc, colour = Father.Acc)) +
theme_bw() +
geom_smooth(method = "gam", formula = y ~ s(x, k = 7, bs = "cs"), se = FALSE) +
geom_ribbon(mapping = aes(x = Challenge, ymin = c(Lower.CI), ymax = c(Upper.CI)), alpha = 0.3, colour = NA) +
scale_fill_manual(values = c("dodgerblue3", "magenta"), labels = c("Cold", "Warm")) +
scale_colour_manual(values = c("dodgerblue3", "magenta"), labels = c("Cold", "Warm")) +
labs(fill = "Father Acclimation \n Temperature", colour = "Father Acclimation \n Temperature") +
theme(axis.text = element_text(size = 8, colour = "black", family = "Noto Sans")) +
labs(x = "Challenge Temperature (°C)", y = expression("Resting MO"[2] * " (mg O"[2] * "/h)")) +
theme(panel.grid.major = element_blank()) +
xlim(c(15, 30))
facet_names <- c(`1` = "Cold Fathers", `2` = "Warm Fathers")
facet_names_new <- c(`1` = "Cold Offspring", `2` = "Warm Offspring")
EM_Plot_Off_Pat_Ribbon <- ggplot(Off_Pat, aes(x = Challenge, y = yvar, fill = Father.Acc, colour = Father.Acc)) +
theme_bw() +
geom_smooth(method = "gam", formula = y ~ s(x, k = 7, bs = "cs"), se = FALSE) +
geom_ribbon(mapping = aes(x = Challenge, ymin = c(Lower.CI), ymax = c(Upper.CI)), alpha = 0.3, colour = NA) +
scale_fill_manual(values = c("dodgerblue3", "magenta"), labels = c("Cold", "Warm")) +
scale_colour_manual(values = c("dodgerblue3", "magenta"), labels = c("Cold", "Warm")) +
labs(colour = "Father Acclimation \n Temperature", fill = "Father Acclimation \n Temperature") +
theme(axis.text = element_text(size = 8, colour = "black", family = "Noto Sans")) +
facet_wrap(~Offspring.Acc, labeller = as_labeller(facet_names_new)) +
labs(x = "Challenge Temperature (°C)", y = expression("Resting MO"[2] * " (mg O"[2] * "/h)")) +
theme(strip.text.x = element_text(family = "Noto Sans", colour = "black"), panel.grid.major = element_blank()) +
xlim(c(15, 30))
facet_names <- c(`1` = "Cold Mothers", `2` = "Warm Mothers")
EM_Plot_Off_Mat_Ribbon <- ggplot(Off_Mat, aes(x = Challenge, y = yvar, fill = Mother.Acc, colour = Mother.Acc)) +
theme_bw() +
geom_smooth(method = "gam", formula = y ~ s(x, k = 7, bs = "cs"), se = FALSE) +
geom_ribbon(mapping = aes(x = Challenge, ymin = c(Lower.CI), ymax = c(Upper.CI)), alpha = 0.3, colour = NA) +
scale_fill_manual(values = c("dodgerblue3", "magenta"), labels = c("Cold", "Warm")) +
scale_colour_manual(values = c("dodgerblue3", "magenta"), labels = c("Cold", "Warm")) +
labs(colour = "Mother Acclimation \n Temperature", fill = "Mother Acclimation \n Temperature") +
theme(axis.text = element_text(size = 8, colour = "black", family = "Noto Sans")) +
facet_wrap(~Offspring.Acc, labeller = as_labeller(facet_names_new)) +
labs(x = "Challenge Temperature (°C)", y = expression("Resting MO"[2] * " (mg O"[2] * "/h)")) +
theme(strip.text.x = element_text(family = "Noto Sans", colour = "black"), panel.grid.major = element_blank()) +
xlim(c(15, 30))
EM_Plot_Mass <- ggplot(EM_Mass, aes(x = Mass, y = yvar)) +
theme_bw() +
geom_smooth(method = "gam", formula = y ~ s(x, k = 7, bs = "cs"), se = FALSE, colour = "black") +
geom_ribbon(mapping = aes(x = Mass, ymin = c(Lower.CI), ymax = c(Upper.CI)), alpha = 0.3, colour = NA) +
xlab("Mass (g)") +
theme(axis.text = element_text(size = 8, colour = "black", family = "Noto Sans")) +
labs(x = "Mass (g)", y = expression("MO"[2] * " (mg O"[2] * "/h)")) +
theme(panel.grid.major = element_blank())
Producing marginal mean plots for glmmTMB model
{
EM_Means_Off <- emmeans::emmip(ModelTMB.AR, Offspring.Acc ~ bSpline, cov.reduce = FALSE)
EM_Means_Mat <- emmeans::emmip(ModelTMB.AR, Mother.Acc ~ bSpline, cov.reduce = FALSE)
EM_Means_Pat <- emmeans::emmip(ModelTMB.AR, Father.Acc ~ bSpline, cov.reduce = FALSE)
EM_Means_Off_Pat <- emmeans::emmip(ModelTMB.AR, Offspring.Acc ~ bSpline | Father.Acc, cov.reduce = FALSE)
EM_Means_Off_Mat <- emmeans::emmip(ModelTMB.AR, Offspring.Acc ~ bSpline | Mother.Acc, cov.reduce = FALSE)
EM_Means_Mass <- emmeans::emmip(ModelTMB.AR, ~Mass, at = list(Mass = seq(min(Dat$Mass, na.rm = T), max(Dat$Mass, na.rm = T), 0.1)), CIs = TRUE, type = "response", data = Dat)
}
Raw <- as.data.frame(emmeans::emmeans(ModelTMB.AR, specs = "bSpline", by = "Offspring.Acc", cov.reduce = FALSE))
bSpline_vals <- Raw$bSpline
Challenge_Temp <- c()
for (i in 1:length(bSpline_vals)) {
Challenge_Temp[i] <- BTChallenge[c(which(BTChallenge$bSpline == bSpline_vals[i])), 9]
}
# Removing replicates again
Challenge_Temp <- unique(Challenge_Temp)
EM_Off <- as.data.frame(EM_Means_Off$data)
EM_Off$Challenge <- sort(rep(Challenge_Temp, 2))
# Removing prodictions for warm acclimated offspring at cold temperatures
to_rm <- c(which(EM_Off$Offspring.Acc == "2" & EM_Off$Challenge < 19))
EM_Off <- EM_Off[-to_rm, ]
EM_Off$Upper.CI <- EM_Off$yvar + EM_Off$SE * 1.96
EM_Off$Lower.CI <- EM_Off$yvar - EM_Off$SE * 1.96
EM_Plot_Off_Ribbon <- ggplot(EM_Off, aes(x = Challenge, y = yvar, fill = Offspring.Acc, linetype = Offspring.Acc)) +
theme_bw() +
geom_smooth(method = "gam", formula = y ~ s(x, k = 7, bs = "cs"), se = FALSE, colour = "black") +
geom_ribbon(mapping = aes(x = Challenge, ymin = c(Lower.CI), ymax = c(Upper.CI)), alpha = 0.3, colour = NA) +
scale_fill_manual(values = c("dodgerblue3", "magenta"), labels = c("Cold", "Warm")) +
scale_linetype_manual(values = c("solid", "dashed"), labels = c("Cold", "Warm")) +
labs(fill = "Offspring Acclimation \n Temperature", linetype = "Offspring Acclimation \n Temperature") +
theme(axis.text = element_text(size = 8, colour = "black", family = "Noto Sans")) +
labs(x = "Challenge Temperature (°C)", y = expression("Resting MO"[2] * " (mg O"[2] * "/h)")) +
theme(panel.grid.major = element_blank()) +
xlim(c(15, 30))
EM_Mat <- as.data.frame(EM_Means_Mat$data)
EM_Mat$Challenge <- sort(rep(Challenge_Temp, 2))
EM_Mat$Upper.CI <- EM_Mat$yvar + EM_Mat$SE * 1.96
EM_Mat$Lower.CI <- EM_Mat$yvar - EM_Mat$SE * 1.96
EM_Plot_Mat_Ribbon <- ggplot(EM_Mat, aes(x = Challenge, y = yvar, linetype = Mother.Acc, fill = Mother.Acc)) +
theme_bw() +
geom_smooth(method = "gam", formula = y ~ s(x, k = 7, bs = "cs"), se = FALSE, colour = "black") +
geom_ribbon(mapping = aes(x = Challenge, ymin = c(Lower.CI), ymax = c(Upper.CI)), alpha = 0.3, colour = NA) +
scale_fill_manual(values = c("dodgerblue3", "magenta"), labels = c("Cold", "Warm")) +
scale_linetype_manual(values = c("solid", "dashed"), labels = c("Cold", "Warm")) +
labs(fill = "Mother Acclimation \n Temperature", linetype = "Mother Acclimation \n Temperature") +
theme(axis.text = element_text(size = 8, colour = "black", family = "Noto Sans")) +
labs(x = "Challenge Temperature (°C)", y = expression("Resting MO"[2] * " (mg O"[2] * "/h)")) +
theme(panel.grid.major = element_blank()) +
xlim(c(15, 30))
EM_Pat <- as.data.frame(EM_Means_Pat$data)
EM_Pat$Challenge <- sort(rep(Challenge_Temp, 2))
EM_Pat$Upper.CI <- EM_Pat$yvar + EM_Pat$SE * 1.96
EM_Pat$Lower.CI <- EM_Pat$yvar - EM_Pat$SE * 1.96
EM_Plot_Pat_Ribbon <- ggplot(EM_Pat, aes(x = Challenge, y = yvar, fill = Father.Acc, linetype = Father.Acc)) +
theme_bw() +
geom_smooth(method = "gam", formula = y ~ s(x, k = 7, bs = "cs"), se = FALSE, colour = "black") +
geom_ribbon(mapping = aes(x = Challenge, ymin = c(Lower.CI), ymax = c(Upper.CI)), alpha = 0.3, colour = NA) +
scale_fill_manual(values = c("dodgerblue3", "magenta"), labels = c("Cold", "Warm")) +
scale_linetype_manual(values = c("solid", "dashed"), labels = c("Cold", "Warm")) +
labs(fill = "Father Acclimation \n Temperature", linetype = "Father Acclimation \n Temperature") +
theme(axis.text = element_text(size = 8, colour = "black", family = "Noto Sans")) +
labs(x = "Challenge Temperature (°C)", y = expression("Resting MO"[2] * " (mg O"[2] * "/h)")) +
theme(panel.grid.major = element_blank()) +
xlim(c(15, 30))
Off_Pat <- as.data.frame(EM_Means_Off_Pat$data)
Off_Pat$Challenge <- sort(rep(Challenge_Temp, 4))
Off_Pat$Upper.CI <- Off_Pat$yvar + Off_Pat$SE * 1.96
Off_Pat$Lower.CI <- Off_Pat$yvar - Off_Pat$SE * 1.96
to_rm <- c(which(Off_Pat$Offspring.Acc == "2" & Off_Pat$Challenge < 19))
Off_Pat <- Off_Pat[-to_rm, ]
facet_names <- c(`1` = "Cold Fathers", `2` = "Warm Fathers")
facet_names_new <- c(`1` = "Cold Offspring", `2` = "Warm Offspring")
EM_Plot_Off_Pat_Ribbon <- ggplot(Off_Pat, aes(x = Challenge, y = yvar, fill = Father.Acc, linetype = Father.Acc)) +
theme_bw() +
geom_smooth(method = "gam", formula = y ~ s(x, k = 7, bs = "cs"), se = FALSE, colour = "black") +
geom_ribbon(mapping = aes(x = Challenge, ymin = c(Lower.CI), ymax = c(Upper.CI)), alpha = 0.3, colour = NA) +
scale_fill_manual(values = c("dodgerblue3", "magenta"), labels = c("Cold", "Warm")) +
scale_linetype_manual(values = c("solid", "dashed"), labels = c("Cold", "Warm")) +
labs(linetype = "Father Acclimation \n Temperature", fill = "Father Acclimation \n Temperature") +
theme(axis.text = element_text(size = 8, colour = "black", family = "Noto Sans")) +
facet_wrap(~Offspring.Acc, labeller = as_labeller(facet_names_new)) +
labs(x = "Challenge Temperature (°C)", y = expression("Resting MO"[2] * " (mg O"[2] * "/h)")) +
theme(strip.text.x = element_text(family = "Noto Sans", colour = "black"), panel.grid.major = element_blank()) +
xlim(c(15, 30))
Off_Mat <- as.data.frame(EM_Means_Off_Mat$data)
Off_Mat$Challenge <- sort(rep(Challenge_Temp, 4))
Off_Mat$Upper.CI <- Off_Mat$yvar + Off_Mat$SE * 1.96
Off_Mat$Lower.CI <- Off_Mat$yvar - Off_Mat$SE * 1.96
to_rm <- c(which(Off_Mat$Offspring.Acc == "2" & Off_Mat$Challenge < 19))
Off_Mat <- Off_Mat[-to_rm, ]
facet_names <- c(`1` = "Cold Mothers", `2` = "Warm Mothers")
EM_Plot_Off_Mat_Ribbon <- ggplot(Off_Mat, aes(x = Challenge, y = yvar, fill = Mother.Acc, linetype = Mother.Acc)) +
theme_bw() +
geom_smooth(method = "gam", formula = y ~ s(x, k = 7, bs = "cs"), se = FALSE, colour = "black") +
geom_ribbon(mapping = aes(x = Challenge, ymin = c(Lower.CI), ymax = c(Upper.CI)), alpha = 0.3, colour = NA) +
scale_fill_manual(values = c("dodgerblue3", "magenta"), labels = c("Cold", "Warm")) +
scale_linetype_manual(values = c("solid", "dashed"), labels = c("Cold", "Warm")) +
labs(linetype = "Mother Acclimation \n Temperature", fill = "Mother Acclimation \n Temperature") +
theme(axis.text = element_text(size = 8, colour = "black", family = "Noto Sans")) +
facet_wrap(~Offspring.Acc, labeller = as_labeller(facet_names_new)) +
labs(x = "Challenge Temperature (°C)", y = expression("Resting MO"[2] * " (mg O"[2] * "/h)")) +
theme(strip.text.x = element_text(family = "Noto Sans", colour = "black"), panel.grid.major = element_blank()) +
xlim(c(15, 30))
EM_Mass <- as.data.frame(EM_Means_Mass$data)
EM_Mass$Upper.CI <- EM_Mass$yvar + EM_Mass$SE * 1.96
EM_Mass$Lower.CI <- EM_Mass$yvar - EM_Mass$SE * 1.96
EM_Plot_Mass <- ggplot(EM_Mass, aes(x = Mass, y = yvar)) +
theme_bw() +
geom_smooth(method = "gam", formula = y ~ s(x, k = 7, bs = "cs"), se = FALSE, colour = "black") +
geom_ribbon(mapping = aes(x = Mass, ymin = c(Lower.CI), ymax = c(Upper.CI)), alpha = 0.3, colour = NA) +
xlab("Mass (g)") +
theme(axis.text = element_text(size = 8, colour = "black", family = "Noto Sans")) +
labs(x = "Mass (g)", y = expression("MO"[2] * " (mg O"[2] * "/h)")) +
theme(panel.grid.major = element_blank())
grid.arrange(EM_Plot_Off_Ribbon, EM_Plot_Mat_Ribbon, EM_Plot_Pat_Ribbon, EM_Plot_Off_Mat_Ribbon, EM_Plot_Off_Pat_Ribbon, EM_Plot_Mass, nrow = 3)
Spline_plot <- BTChallenge %>%
ggplot(aes(x = Challenge.Temp, y = MO2.NMS)) +
geom_point(pch = 21, size = 3, alpha = 0.3, colour = "black", fill = "royalblue3") +
geom_smooth(
mapping = aes(x = Challenge.Temp, y = bSpline),
colour = "black", linetype = "dashed"
) +
labs(x = "Challenge Temperature (°C)", y = expression("MO"[2] * " (mg O"[2] * "/h)")) +
theme_bw() +
theme(axis.text = element_text(size = 8, colour = "black", family = "Noto Sans"))
print(Spline_plot)
dev.new()
EM_Plot_Off_Ribbon
showtext_opts(dpi = 800)
showtext_auto()
ggsave("/home/joshk/Documents/Trout/Brook Trout Challenge Data/Offspring_alone_z.jpeg", EM_Plot_Off_Ribbon, dpi = 800, width = 9.45, height = 6.72)
showtext_auto(enable = FALSE)
dev.new()
EM_Plot_Mat_Ribbon
showtext_opts(dpi = 800)
showtext_auto()
ggsave("/home/joshk/Documents/Trout/Brook Trout Challenge Data/Mothers_alone_z.jpeg", EM_Plot_Mat_Ribbon, dpi = 800, width = 9.45, height = 6.72)
showtext_auto(enable = FALSE)
dev.new()
EM_Plot_Pat_Ribbon
showtext_opts(dpi = 800)
showtext_auto()
ggsave("/home/joshk/Documents/Trout/Brook Trout Challenge Data/Fathers_alone_z.jpeg", EM_Plot_Pat_Ribbon, dpi = 800, width = 9.45, height = 6.72)
showtext_auto(enable = FALSE)
dev.new()
EM_Plot_Off_Mat_Ribbon
showtext_opts(dpi = 800)
showtext_auto()
ggsave("/home/joshk/Documents/Trout/Brook Trout Challenge Data/Offspring_by_mothers_z.jpeg", EM_Plot_Off_Mat_Ribbon, dpi = 800, width = 10, height = 6.72)
showtext_auto(enable = FALSE)
dev.new()
EM_Plot_Off_Pat_Ribbon
showtext_opts(dpi = 800)
showtext_auto()
ggsave("/home/joshk/Documents/Trout/Brook Trout Challenge Data/Offspring_by_fathers_z.jpeg", EM_Plot_Off_Pat_Ribbon, dpi = 800, width = 10, height = 6.72)
showtext_auto(enable = FALSE)
dev.new()
EM_Plot_Mass
showtext_opts(dpi = 800)
showtext_auto()
ggsave("/home/joshk/Documents/Trout/Brook Trout Challenge Data/Mass_alone_z.jpeg", EM_Plot_Mass, dpi = 800, width = 9.45, height = 6.72)
showtext_auto(enable = FALSE)
dev.new()
Spline_plot
showtext_opts(dpi = 800)
showtext_auto()
ggsave("/home/joshk/Documents/Trout/Brook Trout Challenge Data/Spline_z.jpeg", Spline_plot, dpi = 800, width = 9.45, height = 6.72)
showtext_auto(enable = FALSE)